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Abstract 


The  research  addressed  interactions  amongst  droplets  in  a  dense  spray.  The  effects 
of  neighboring  droplets,  that  were  a  few  droplet  diameters  away,  on  a  vaporizing  droplet 
were  examined  by  theoretical  and  computationcil  analyses  for  two  basic  configurations: 
(1)  the  axisymmetric  convective  situation  where  two  or  three  droplets  moved  in  tandem 
and  (2)  the  fully  three-dimensional  convective  situation  where  droplets  moved  side-by- 
side.  Droplets  in  the  wake  of  other  droplets  experienced  a  reduction  in  drag  force,  trans¬ 
port  rates,  and  vaporization  rate,  sometimes  causing  collisions.  Sufficiently  close  droplets 
moving  side-by-side,  approximately  in  parallel,  experienced  a  repulsive  lift  force  and  an 
increased  drag  force.  Vaporizing  liquid  oxygen  droplets  in  a  hydrogen  gas  environment 
were  studied  at  both  subcritical  and  supercritical  pressures  considering  the  variable  liquid 
density  with  the  associated  droplet  swelling  during  heating  and  the  dependence  of  the  lo¬ 
cal  critical  state  upon  local  composition.  Droplet  surface  conditions  could  be  subcritical 
even  if  pressures  were  supercritical  for  pure  oxygen  due  to  diffusing  hydrogen.  The  critical 
surface  regressed  towards  the  droplet  surface  as  the  droplet  heated.  Engineering  correla¬ 
tions  for  the  drag  coefficients,  Nusselt  numbers,  and  Sherwood  numbers  for  hydrocarbon 
fuel  droplets  in  dense  sprays  were  obtained. 
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N  o  menclat  ure 


Co  initial  droplet  radius 

at  instantaneous  droplet  radius 

do  initial  distance  between  droplet  centers  normalized  initial  droplet  diameter 

dt  instantaneous  distance  between  droplet  centers  normalized  initial  droplet  diameter 

D'oo  dimensional  mass  difFusivity  at  infinity 

Dij  binary  diffusion  coefficient 

/  fugacity 

h'  dimensional  film  heat  transfer  coefficient 

h'^  dimensional  film  mass  transfer  coefficient 

AHy  enthalpy  of  vaporization 

k'^  dimensional  thermal  conductivity  at  infinity 

k  thermal  conductivity  of  the  gas  phcise 

rhy  mass  vaporization  rate  on  the  droplet  surface 

yf,av  average  mass  fraction  on  the  droplet  surface 

i,  j,  k  unit  vectors  in  x,  y,  and  z  direction,  respectively 
It  instantaneous  droplet  position  in  y  direction 

numbers  of  grid  points  in  directions 

dimensional  kinematic  viscosity  at  infinity 
P  pressure 

r  radial  direction 

Re  initial  Reynolds  number  based  on  initial  droplet  diameter,  U^o^a'^f 

Ret  instantaneous  Reynolds  number, 

T  temperature 

t  dimensionless  time  normalized  by 

U'^  g  initial  dimensional  droplet  velocity 

U'^  t  instantaneous  dimensional  droplet  velocity,  yjv^t  + 

U'gg  t  instantaneous  dimensional  droplet  velocity  in  x-direction 

V^,t  instantaneous  dimensional  droplet  velocity  in  y-direction 

X,  y,  z  Cartesian  coordinates 

X  mole  fraction 

Y  mass  fraction 


Greek  symbols 

p  dimensionless  gas  density  normalized  by 

p'gg  dimensional  gas  density  at  infinity 

6  angular  direction 

T},  (  nonorthogonal  generalized  coordinates 
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■o-  -e-  H 


viscous  stress  tensor 
pseudo- stream  function 
fugacity  coefficient 


Superscript 

'  dimensional  quantity 


Subscript 

o  initial  qucintity 

/  liquid  phase 


1.  Introduction 


Practical  combustors  with  direct  injection  of  liquid  fuel  have  regions  of  large  con¬ 
centration  of  droplets,  which  we  will  call  dense  sprays.  In  such  regions,  the  effect  of 
neighboring  droplets  modifies  the  ambient  conditions  in  the  flow  near  any  given  droplet. 
As  a  consequence,  the  drag  coefficient,  lift  coefficient,  moment  coefficient,  Nusselt  number, 
Sherwood  number,  and  vaporization  rates  are  different  from  those  of  an  isolated  droplet 
at  the  same  Reynolds  number  and  transfer  number.  Also,  the  droplet  trajectories,  droplet 
lifetime,  and  local  gas-phase  mixture  ratios  can  be  significantly  affected  by  the  differences 
in  transport  processes  and  the  modified  flow  field  due  to  neighboring  droplets.  In  many 
practical  situations,  including  hquid  rocket  engines,  turbojet  engines,  and  ramjet  engines, 
droplet  heating  and  vaporization  can  be  a  rate-controlling  phenomena.  Therefore,  the 
modifications  in  the  droplet  behavior  due  to  dense  spray  effects  can  have  important  con¬ 
sequences. 

The  geometrical  configurations  of  droplets  in  a  real  region  of  large  concentration  are 
complex  and  subject  to  uncertainty.  Droplet  arrays  as  discussed  by  Sirignano  [1],  although 
artificial,  can  provide  information  on  droplets  interaction  and  give  a  detailed  analysis  of 
the  problem.  The  principal  investigator  and  his  co-workers  [2-7]  performed  research  on 
axisymmetric  configurations  of  interacting  particles  (or  droplets)  in  the  wake  of  another 
particles  (or  droplets).  The  major  new  goal  of  this  project  is  the  study  of  three-dimensional 
interactions  between  droplets. 

Three  tasks  hcis  been  performed  step  by  step  to  attack  the  three-dimensional  in¬ 
teractions  between  droplets.  First,  we  investigated  incompressible,  isothermal,  three- 
dimensional  flow  over  two  identical  (solid  or  liquid)  spheres  which  were  held  fixed  relative 
to  each  other  in  the  transverse  direction  against  the  uniform  stream  at  Reynolds  number 
0(100).  We  determined  the  effects  of  three-dimensional  interactions  on  the  lift,  moment, 
and  drag  coefficients  as  a  function  of  the  dimensionless  distance  between  the  two  spheres 
and  Reynolds  number.  Some  novel  phenomena  in  the  near  wake  were  discovered  as  the 
gap  between  the  two  spheres  decreases.  Our  results  are  also  compared  with  available  ex¬ 
perimental  and  numerical  data.  This  study  will  be  discussed  briefly  in  Section  2.  Details 
are  provided  in  Reference  9  which  is  attached  as  an  addendum  to  this  report. 

Secondly,  we  studied  interactions  of  two  identical  droplets  which  are  injected  and  then 
move  side  by  side  into  an  initially  quiescent,  incompressible,  isothermal  fluid  medium. 
The  droplets  decelerate  due  to  the  drag  and  change  their  direction  of  motion  due  to  the 
lift.  By  placing  the  origin  of  a  noninertial  reference  frame  at  the  center  of  mass  of  the 
two  droplets,  the  Navier-Stokes  equations  include  a  noninertial  term  which  is  evaluated 
from  Newton’s  second  law  for  the  droplet  motion.  This  study  will  be  discussed  briefly  in 
Section  3.  Details  are  provided  in  Reference  18  which  is  attached  as  an  addendum  to  this 
report. 

In  the  third  task,  we  extended  the  above  studies  by  considering  three-dimensional 
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interactions  between  two  identical  n-octane  droplets  injected  into  a  high-temperature  and 
high-pressure  environment  which  is  typically  encountered  in  practical  combustors.  The 
droplets  are  cold  initially  and  are  subsequently  vaporizing  due  to  the  high  temperature. 
The  energy  equation,  the  species  equation,  and  the  equation  of  state,  in  addition  to 
the  Navier-Stokes  equations,  are  solved  simultaneously  in  the  gas  phase  to  account  for 
vaporization  of  the  droplets  and  temperature  distribution  in  the  domain.  In  addition, 
nonlinear  coupled  interface  boundary  conditions  between  the  gas  and  liquid  phases  are 
imposed  simultaneously.  This  study  will  be  discussed  in  detail  in  Section  4.  We  expect 
to  submit  a  manuscript  for  journal  publication  on  this  subject  in  the  near  future. 

In  addition  to  the  above  three-dimensional  calculations,  three  tasks  which  make  use 
of  advanced  axisymmetrical  models  have  been  conducted  in  this  program,  namely:  (1) 
axisymmetric  interactive  fuel  droplet  calculations,  (2)  LOX  vaporization  at  subcritical 
conditions,  (3)  LOX  vaporization  at  near-critical  or  supercritical  conditions. 

The  first  task  investigates  fuel  droplets  moving  in  tandem  by  axisymmetric  calcu¬ 
lations.  This  research  addresses  the  interaction  of  three  vaporizing  droplets  moving 
collinearly  which  represents  a  model  of  an  injected  stream  of  fuel  droplets.  The  pur¬ 
poses  of  this  study  are  to  study  the  wake  effect  of  the  lead  droplet  on  the  downstream 
droplets  and  to  examine  the  effects  of  initial  spacing  on  the  total  system.  Some  technical 
discussion  will  be  given  in  Section  5.  More  details  are  available  in  our  journal  publication 
[7].  This  publication  is  attached  as  an  addendum  to  this  report. 

In  another  task,  vaporizing  liquid  oxygen  droplets  are  studied  including  axisymmetric 
interactions  amongst  fuel  and  oxygen  droplets.  Vaporization  of  LOX  droplets  in  convective 
hydrogen  environment  over  a  wide  range  of  pressure  has  been  investigated.  Current 
emphasis  is  placed  on  the  understanding  of  supercritical  LOX  vaporization  where  many 
complex  transport  phenomena,  such  cis  the  vanishing  of  the  gas/liquid  interface,  diffusion 
of  gas  vapor  into  droplet,  and  real  gas  effects,  become  very  important  in  determining  the 
droplet  vaporization  rate. 

Progress  has  occurred  in  the  study  of  vaporization  of  an  isolated  LOX  droplet  in  a 
convective  hydrogen  environment  at  subcritical  conditions.  This  study  will  be  discussed 
in  Section  6.  A  copy  of  the  preprint  describing  this  research  is  attached  as  an  addendum 
to  this  report.  Modifications  of  this  low  pressure  model  to  consider  near-critical  or  super¬ 
critical  conditions  have  been  made  in  order  to  fulfill  this  task.  The  details  are  presented 
in  Section  7. 


2.  Three-dimensional  flow  over  two  spheres  placed  side  by  side 

2.1  Problem  statement  andi  numerical  solution 

We  consider  a  steady,  three-dimensional,  incompressible,  laminar  flow  of  a  Newtonian 
fluid  past  two  identical  (solid  or  liquid)  spheres  held  fixed,  with  the  line  connecliug  the 


5 


sphere  centers  normal  to  a  uniform  stream,  as  shown  in  Figure  1;  do  denotes  the  distance, 
normalized  by  the  sphere  radius,  from  the  sphere  center  to  the  x-y  symmetry  plane  be¬ 
tween  the  two  spheres.  Far  upstream,  the  flow  is  uniform  with  constant  velocity  Uooi 
parallel  to  the  x-axis.  Two  symmetry  planes  are  noted  in  Figure  1:  the  x-z  plane  contain¬ 
ing  the  centers  of  the  two  spheres  and  the  x-y  plane  located  at  z  =  —do  midway  between 
the  sphere  centers. 

The  governing  equations  and  boundary  conditions,  coordinate  system,  computational 
grid,  and  numerical  method  are  discussed  in  detail  in  the  addendum  to  this  report. 

The  dimensional  drag,  lift,  and  moment  coefficients  are  evaluated  from 


II 

1  -p'n  -i  dS'  +  f  n-r'-i  dS\ 

's  Js 

(1) 

^  — p'n*k  dS*  +  [  n*r'*k  dS\ 

s  Js 

(2) 

=  j^r'x  r'd5'. 

(3) 

where  S*  denotes  the  surface  of  the  sphere,  n  is  the  outward  unit  normal  vector  at  the 
surface,  r'  is  the  position  vector  from  the  center  of  the  sphere,  and  r'  is  the  viscous 
stress  tensor.  The  lift  force  is  assumed  positive  when  it  is  directed  towaxd  the  positive 
z-axis.  Due  to  symmetry,  only  the  y-component  of  the  moment  is  non-zero  and  is  assumed 
positive  in  counter-clockwise  direction. 

Nondimensional  coefficients  of  drag,  lift,  and  moment  are  defined  respectively  as 


Cd  = 

F'd 

(4) 

Cl  = 

n 

(5) 

Cm  = 

M'-j 

(6) 

2.2  Results  and  discussion 

We  first  tested  the  three-dimensional  code  by  solving  the  flow  generated  by  an  impul¬ 
sively  started  solid  sphere  in  a  quiescent  fluid  at  two  Reynolds  numbers:  50  and  100.  The 
time-dependent  solution  converges  asymptotically  to  a  steady-state  which  is  in  excellent 
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agreement  with  available  experimental  data  and  correlations  as  shown  in  Tables  1,  2,  3, 
and  4.  Tables  1  and  3  list  the  drag  coefficients  as  a  function  of  the  computational  grid 
density  at  Reynolds  numbers  50  and  100  respectively,  and  compare  them  with  the  data 
of  Roos  &  Willmarth  [10]  and  also  with  the  correlations  of  Clift  et  al.  [11].  Tables  2  and  4 
show  the  pressures  at  the  front  and  rear  stagnation  points  and  the  separation  angle  mea¬ 
sured  from  the  front  stagnation  point  as  a  function  of  grid  density  at  Reynolds  number  50 
and  100,  in  comparison  with  the  data  of  Taneda  [12]  and  also  with  the  correlations  of  Clift 
et  al.  [11].  Although  the  solution  in  these  test  cases  are  axisymmetric,  none  of  the  three 
velocity  components  in  our  formulation  becomes  identically  zero.  Therefore,  the  three- 
dimensional  solution  scheme  is  fully  exercised  here.  The  calculations  were  performed  for 
three  different  grids,  {Ni  x  N^y.  N^)  =  (20  x  21  x  21),  (30  x  31  x  31),  and  (40  x  41  x  41), 
in  a  computational  domain  with  an  outer  boundary  located  at  21  sphere  radii  from  the 
sphere  center.  The  dimensionless  times  needed  to  reach  steady  state  for  Reynolds  number 

50  and  100  are  9  and  15,  respectively.  The  results  for  interactions  of  two  solid  spheres  will 
be  presented  and  discussed  first  and,  then,  those  for  two  hquid  spheres  will  be  discussed 
shortly  thereafter  in  the  remainder  of  the  section. 

In  order  to  illustrate  better  the  fluid  motion,  we  consider  the  flow  field  in  the  x-z  plane 
of  symmetry,  which  is  defined  as  the  principal  plane,  where  the  narrowest  path  between 
the  two  spheres  occurs;  hence  the  strongest  interactions  between  them  occur  there.  We 
first  describe  the  main  flow  features  that  emerged  from  our  results,  then  we  examine  the 
results  in  detail. 

Figure  2(a)  displays  a  sketch  of  actual  streamlines  around  a  single  solid'  or  hquid 
sphere  in  the  principal  plane  at  Reynolds  number  100.  As  expected,  two  identical  counter¬ 
rotating  vortices  (corresponding  to  a  single  vortex  ring)  exist  in  the  wake,  and  the  down¬ 
stream  stagnation  point  is  located  on  the  axis  of  symmetry.  Line  I-II  connecting  the  front 
and  rear  stagnation  points  in  the  standard  axisymmetric  flow  over  a  single  sphere  will  be 
used  as  a  reference  line.  We  refer  to  the  region  above  the  line  I-II  as  the  ‘top’  or  ‘upper’ 
region  and  that  below  the  Une  as  the  ‘bottom’  or  ‘lower’  region. 

Figure  2(b)  displays  a  sketch  of  actual  streamhnes  around  one  of  the  two  spheres  in 
the  principal  plane  at  Reynolds  number  100.  The  two  spheres  are  separated  by  a  distance 
do  —  2.  Due  to  the  blockage  of  the  flow  in  the  gap  between  the  two  spheres,  the  streamlines 
diverge  away  from  the  x-y  symmetry  plane  (located  at  ^  =  —do)  cis  they  approach  the 
front  stagnation  region.  Thus,  the  stagnation  streamline  of  the  single  sphere  case  (I- 

51  in  Figure  2(a))  no  longer  intersects  the  sphere,  and  another  streamline  closer  to  the 
symmetry  plane  meets  the  sphere  to  form  the  new  front  stagnation  streamline  at  point 
P.  As  a  consequence,  the  fluid  particles  move  faster  in  the  lower  left  region  around  the 
sphere  than  in  the  upper  left  region,  and  this  causes  the  pressure  in  the  lower  left  region 
to  be  lower  than  that  in  the  upper  left  region.  The  resulting  pressure  difference  between 
the  upper  and  lower  left  regions  is  higher  than  that  between  the  bottom  of  the  sphere  and 
the  narrow  path.  This  pressure  imbalance,  which  will  be  discussed  later  in  detail,  causes 
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repulsion  of  the  two  spheres.  The  contribution  of  shear  stress  differences  to  the  repulsion 
will  be  discussed  later  in  this  section. 

Figure  2(b)  shows  an  interesting  streamline  pattern  in  the  wake  region.  Two  counter¬ 
rotating  eddies  exist  in  the  wake  but  their  configuration  is  quite  different  from  that  for 
axisymmetric  flow.  The  lower  eddy  is  formed  by  the  fluid  separating  on  the  lower  portion 
of  the  sphere  as  in  the  case  of  axisymmetric  flow.  The  upper  eddy  is  not  formed  by 
the  fluid  separating  on  the  upper  portion  of  the  sphere,  but  rather  by  the  fluid  turning 
around  the  lower  eddy  and  being  entrained  by  the  upper  flow.  This  upper  eddy  is  detached 
from  the  sphere.  A  portion  of  the  fluid  moving  around  the  bottom  of  the  sphere  passes 
between  the  detached  upper  eddy  and  the  sphere.  The  streamline  A-B  encompassing  the 
upper  eddy  intersects  itself,  and  the  intersection  point,  C,  designated  as  the  downstream 
stagnation  point,  is  shifted  toward  the  x-y  symmetry  plane.  Both  eddies  are  smaller  than 
those  of  the  axisymmetric  flow.  These  new  features  can  be  explained  as  follows.  The 
pressure  above  the  upper  wake  is  less  than  that  below  the  lower  wake  due  to  the  increased 
acceleration  of  the  fluid  in  the  narrow  path  between  the  two  spheres.  Thus,  the  fluid 
particles  turning  around  the  lower  eddy  are  pushed  into  the  upper  region  of  the  wake. 
The  pressure  distribution  around  the  sphere  will  be  discussed  further. 

Figure  2(c)  shows  a  sketch  of  the  actual  streamline  pattern  at  Reynolds  number  100 
for  the  case  of  do  =  1.5.  The  shifting  of  the  front  stagnation  streamline  and  stagnation 
point  toward  the  x-y  symmetry  plane  are  more  obvious  here  than  in  the  previous  case 
of  do  =  2.  The  significant  difference  is  in  the  wake  region  where  both  the  upper  eddy 
and  downstream  stagnation  point  vanish.  Fluid  particles  separating  on  the  uppfer  portion 
of  the  sphere  move  downstream  without  returning  (streamline  D-E).  On  the  other  hand, 
fluid  particles  turning  around  the  lower  eddy  move  into  the  upper  region  of  the  wake  until 
they  reach  near  the  upper  separation  point,  D,  and  then  move  downstream  in  an  S-shaped 
path  (streamline  F-G)  without  returning  to  form  an  eddy.  The  lower  eddy  shrinks  as  the 
two  spheres  are  closer,  and  the  pressure  difference  between  the  top  and  bottom  of  the 
wake  is  larger. 

It  is  interesting  to  examine  the  changes  in  the  separation  region  at  the  sphere  surface 
for  the  cases  of  do  =  2  and  1.5.  More  specifically,  we  examine  the  behavior  of  the  circle 
of  intersection  of  the  wake  and  the  sphere.  Our  results  show  that  the  circle  is  slightly 
shifted  toward  the  x-y  symmetry  plane  due  to  the  decreased  pressure  in  the  gap  region 
with  respect  to  that  in  the  wake  lower  region.  This  shifting  produces  separation  angles 
at  the  top,  middle,  and  bottom  of  the  sphere  with  values  of  123.1°,  126.5°,  and  126.2°, 
respectively,  for  the  case  of  do  =  1.5,  where  the  angles  are  measured  from  the  front 
stagnation  point  of  the  axisymmetric  flow  case.  This  is  in  contrast  with  an  angle  of 
125.7°  at  all  positions  in  the  single  sphere  case. 

A  stream  function  cannot  be  defined  and  calculated  from  the  velocity  in  the  principal 
plane  due  to  the  existence  of  a  divergence  associated  with  the  third  component  of  velocity. 
Nevertheless,  for  descriptive  purposes  only,  it  is  convenient  to  use  the  algorithm 
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to  present  approximations  to  the  streamline  pattern.  It  is  understood  that  since  a  true 
streamfunction  does  not  exist,  the  pseudo-streamfunction,  is  dependent  upon  the 
integration  path.  The  above  algorithm  specifically  involves  only  radial  integration;  ug 
can  be  recovered  by  differentiation  of  this  function,  but  Ur  cannot  be  recovered.  The 
streamlines  presented  in  Figures  2(a),  2(b),  and  2(c)  were  based  on  this  algorithm. 

Phenomena  in  the  wake  similar  to  those  described  above  have  been  found  in  a  few 
previous  studies.  Rosfjord  [13]  obtained  results  similar  to  those  in  Figures  2(b)  and  2(c) 
in  his  experimental  and  numerical  studies  of  the  recirculating  flow  region  between  two- 
dimensional,  parallel,  separated  jets.  He  found  that  for  velocity  ratios  between  two  jets 
equal  to  1.11  and  1.25  (upper  to  lower),  two  eddies  exist  near  the  injector  face,  but  the 
upper  eddy  is  detached  from  the  injector  face,  and  a  portion  of  the  fluid  originating  at 
the  lower  jet  is  entrained  by  the  upper  jet,  passing  between  the  detached  upper  eddy  and 
the  injector  face.  He  also  reported  that  for  a  velocity  ratio  of  1.4,  only  the  lower  eddy 
existed  and  a  complete  entrainment  of  the  weaker  jet  was  observed.  In  particular,  the 
flow  S-shaped  loop  near  the  stronger  jet  was  clearly  indicated.  Recently,  Dandy  Dwyer 
[14]  also  found  a  similar  flow  pattern  in  their  numerical  study  of  steady,  uniform  shear 
flow  past  a  single  solid  sphere. 

In  order  to  facilitate  the  visualization  of  the  three-dimensional  character  of  the  flow 
in  the  wake  region  discussed  above,  we  present  the  pathlines  of  selected  fluid  particles  in 
the  addendum. 

In  the  following  discussion,  we  classify  the  proximity  of  the  two  spheres  into  three 
regimes;  close,  intermediate,  and  far-separated,  depending  on  the  values  of  do  and  Re. 
Figures  3  and  4  pertain  to  flows  around  solid  spheres,  while  Figures  5,  6,  7,  and  8  relate 
to  flows  around  hquid  spheres. 

Figures  3(a),  3(b),  and  3(c)  show  the  total  lift  coefficient  and  the  lift  coefficients  due  to 
viscous  and  pressure  contributions,  respectively,  as  a  function  of  do  at  Reynolds  numbers 
50,  100,  and  150.  The  total  lift  coefficient.  Figure  3(a),  is  positive  when  the  two  spheres 
are  close  (  do  <  7.9  for  Re  =  50,  <  4  for  Re  =  100,  and  do  <  3.4  for  Re  =  150  ).  That 

is,  the  two  spheres  repel  each  other,  and  the  repulsion  is  stronger  the  closer  they  are.  Our 
results  show  that  both  the  viscous  and  pressure  contributions  have  an  important  effect 
on  the  repelling  force,  but  the  pressure  contribution  is  more  dominant  when  Re  >  100 
(compare  Figures  3(b)  and  3(c)).  On  the  other  hand,  the  total  lift  coefficient  is  negative 
and  relatively  small  -  that  is,  the  two  spheres  attract  each  other  weakly  -  at  intermediate 
separation  distances  (  7.9  <  do  <  21  for  Re  ~  50,  4  <  do  <  21  for  Re  =  100,  and 
3.4  <  do  <  21  for  Re  =  150  ).  At  these  distances,  the  pressure  is  the  main  contributor 
to  the  attraction  force  at  all  Reynolds  numbers.  The  smaller  the  Reynolds  number,  the 
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smaller  the  pressure  effect,  the  weaker  the  attraction,  and  the  narrower  the  region  of 
attraction  is.  When  do  >  21,  however,  the  lift  vanishes,  and  the  two  spheres  have  no 
interactions  at  any  Reynolds  numbers. 

As  discussed  earlier,  when  the  two  spheres  are  in  close  proximity,  the  front  stagnation 
point  is  shifted  toward  the  x-y  symmetry  plane,  and  the  fluid  particles  accelerate  faster 
in  the  lower  left  region  than  they  accelerate  in  the  upper  left  region  of  the  sphere.  This 
difference  in  acceleration  results  in  a  net  negative  pressure  gradient  normal  to  and  away 
from  the  x-y  symmetry  plane,  contributing  to  the  repulsion  between  the  two  spheres. 
The  shear  stress  is  also  higher  in  the  lower  left  region  than  in  the  upper  left  region. 
Furthermore,  the  shear  force  in  the  lower  left  region,  due  to  its  inclination  with  the  x- 
axis,  contributes  to  both  the  hft  (parallel  to  the  z-axis)  and  drag  (parallel  to  the  x-axis), 
whereas  the  shear  force  at  the  top  of  the  sphere  contributes  mainly  to  the  drag.  Therefore, 
both  the  pressure  and  shear  forces  contribute  to  a  positive  lift  force  (i.e.  the  two  spheres 
repel  each  other)  when  the  two  spheres  are  close. 

On  the  other  hand,  when  the  two  spheres  are  in  the  intermediate  separation  regime,  the 
velocity  vector  distributions  show  that  the  front  stagnation  streamline  is  almost  straight, 
and  thus  the  flow  in  the  lower  left  is  not  affected  by  the  presence  of  the  other  sphere. 
Nevertheless,  the  gap  between  the  two  spheres  causes  the  flow  to  accelerate  slightly  faster 
on  the  top  of  the  sphere  than  on  the  bottom,  and  as  a  result,  the  average  pressure  in  the 
lower  region  is  slightly  higher  than  that  in  the  narrow  gap.  Thus,  the  two  spheres  attract 
each  other  weakly,  and  the  attraction  force  is  mainly  due  to  the  pressure  distribution. 
The  shear  force,  nearly  parallel  to  the  x-axis  at  the  top  of  the  sphere,  contributes  mainly 
to  the  drag  but  not  to  the  lift. 

Figure  4  shows  the  drag  coefficient  as  a  function  of  the  dimensionless  distance  at 
Reynolds  numbers  50,  100,  and  150.  The  drag  increases  with  decreasing  do  when  do  is  less 
than  4  for  all  Reynolds  numbers.  It  increases  slightly  with  increasing  do,  at  intermediate 
separation  distances,  and  eventually  tends  to  that  of  a  single  sphere  when  do  >  21.  The 
drag  increases  as  the  two  spheres  get  close  because  the  shear  stress  on  the  sphere  is 
increased  and  the  pressure  distribution  is  changed  due  to  the  flow  acceleration  on  the 
lower  left  region  as  well  as  in  the  gap  between  them.  Results  for  the  moment  coefficient 
are  provided  in  the  addendum. 

In  the  analysis  of  the  flow  field  past  two  liquid  spheres,  we  use  a  viscosity  ratio  (in¬ 
ternal  to  external  fluid)  of  25  and  density  ratio  of  300.  These  values  are  typical  of  a 
liquid-hydrocarbon  fuel  in  a  high- temperature,  high-pressure  surrounding  gas  generally 
encountered  in  gas  turbine  combustors  (Raju  &  Sirignano  [5]). 

As  in  the  solid  sphere  case,  we  examine  the  flow  field  for  two  liquid  spheres  in  the  x-z 
plane  of  symmetry,  the  principal  plane,  where  the  narrowest  path  between  liquid  spheres 
is  encountered.  The  streamlines  discussed  earlier  in  Figures  2(a),  2(b),  and  2(c)  can  also 
represent  typical  streamlines  in  the  external  flow  of  liquid  spheres.  However,  there  are 
differences  from  the  solid  sphere  case.  First,  the  angle,  measured  from  <?  =  0,  at  which 
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sepaxation  occurs  on  the  sphere  surface  is  much  higher  than  that  of  the  solid  sphere. 
Secondly,  a  closer  examination  of  the  velocity  vector  plot  (Figures  5(a)  and  6(a))  in 
the  wake  region  indicates  that  the  separating  streamline,  instead  of  being  nearly  normal 
to  the  sphere  surface,  now  curves  closer  to  the  sphere  surface  producing  a  “squashed” 
recirculation  zone.  The  length  of  the  recirculating  eddy  is  also  slightly  smaller  than  that 
of  the  solid  sphere. 

Figures  5(a)  and  5(b)  show  the  velocity  vector  fields  of  the  external  and  internal  flows, 
respectively,  in  the  principal  plane  at  Reynolds  number  100  where  the  two  spheres  are 
separated  by  a  distance  do  =  2.  A  secondary  eddy  in  the  liquid-sphere  stern  region  is 
evident  in  both  the  upper  and  lower  regions  in  the  principal  plane,  but  the  eddy  centers 
in  both  regions  are  asymmetrical  with  respect  to  the  z  =  0  plane.  Also,  these  eddies 
are  concomitant  with  the  occurrence  of  the  eddies  in  both  regions  in  the  external  flow. 
Figures  6(a)  and  6(b)  show  the  velocity  vector  fields  of  the  external  and  internal  flows 
in  the  principal  plane  at  Reynolds  number  100  for  the  case  of  do  =  1.5.  The  secondary 
internal  eddy  in  the  liquid-sphere  stern  region  exists  only  in  the  lower  region,  and  the 
secondary  eddy  in  the  upper  region  no  longer  exists.  The  vanishing  secondary  internal 
eddy  in  the  upper  region  is  accompanied  by  the  disappearance  of  the  recirculating  eddy 
in  the  upper  region  in  the  external  flow. 

Calculations  of  the  lift,  moment,  and  drag  coefficients  were  performed  for  dimensionless 
distances  from  the  liquid  sphere  center  to  the  symmetry  plane  between  two  liquid  spheres 
in  the  range  1 .5  <  do  <  25,  for  a  viscosity  ratio  (internal  to  external  fluid)  of  25  and  density 
ratio  of  300  at  Reynolds  numbers  50,  100,  and  150.  Figures  7  and  8  show  the  Coefficients 
of  total  lift  and  drag  as  a  function  of  dimensionless  distance  at  Reynolds  numbers  50,  100, 
and  150.  The  coefficients  of  total  lift,  moment,  and  drag  are  slightly  smaller  in  absolute 
value  than  those  of  the  solid  spheres  at  both  the  repelling  and  attraction  separation 
distances  and  at  all  Reynolds  numbers.  The  lower  coefficients  of  the  liquid  sphere  are 
attributed  to  the  surface  motion  of  the  liquid  sphere  which  reduces  the  velocity  gradient 
and  friction  force.  A  smaller  drag  coefficient  for  the  liquid  sphere  in  axisymmetric  flow 
has  been  also  found  in  earlier  calculations  (Clift  et  al.  [11]). 

We  did  not  perform  calculations  for  Reynolds  number  higher  than  150,  because  it  is 
known  that  the  wake  becomes  unstable  at  a  Reynolds  number  in  the  range  of  130-220 
for  a  single  sphere  (Taneda  [12];  Goldburg  &  Florsheim  [15];  Willmarth  Sz  Roos  [10]; 
Nakamura  [16];  Kim  &  Pearlstein  [17])  and  our  present  goal  is  to  obtain  steady-state 
solutions.  For  the  solution  of  a  flow  including  the  three-dimensional  unsteady  wake,  a 
complete  computational  domain  (i.e.  encompassing  the  two  spheres  without  a  symmetry 
plane)  and  periodic  boundary  conditions  in  ^-direction  will  be  necessary. 
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3.  Three-dimensional  flow  computation  for  two  interacting, 
moving  droplets 

3.1  Problem  statement 

We  consider  an  unsteady,  three-dimensional,  incompressible,  laminar  flow  generated 
by  two  identical  (spherical)  droplets  injected  into  initially  quiescent  constant  property 
Newtonian  fluid.  The  droplets  move  side  by  side  in  the  same  plane  as  shown  in  Figure  9, 
where  dt  denotes  the  instantaneous  distance,  normalized  by  the  droplet  radius,  from  the 
droplet  center  to  the  y-z  symmetry  plane  between  the  two  droplets.  We  neglect  the  net 
gravity  force  acting  on  the  droplet  and  also  assume  that  the  Weber  number  is  small  enough 
that  the  droplet  remains  spherical.  We  choose  the  origin  of  a  non-rotating  noninertial 
reference  frame  at  the  center  of  mass  of  the  two  droplets  and  the  angle  of  initial  droplet 
motion  Oq  =  0.  Far  upstream,  the  flow  is  uniform  and  has  an  instantaneous  velocity  Voo,t  j 
parallel  to  the  y-ajcis.  It  is  noted  that  there  are  two  S3rmmetry  planes  in  Figure  9.  One  is 
the  x-y  plane  in  which  the  centers  of  the  two  droplets  lie,  and  the  other  is  the  y-z  plane 
which  is  midway  between  the  droplet  centers. 

As  shown  in  Figure  9,  we  utilize  three  coordinate  systems  in  our  formulation,  the 
Cartesian  coordinates  (x,y,z)  for  the  gas  phase  whose  origin  is  at  the  center  of  mass  of 
the  two  droplets,  the  Cartesian  coordinates  (x;,  t/i,  zj)  for  the  liquid  phase  whose  origin  is 
at  the  droplet  center,  and  also  the  nonorthogonal  generalized  coordinates  i  is 

the  radial,  t\  is  the  angular,  and  ^  is  the  azimuthal  direction.  Due  to  symmetry  of  the 
geometry,  the  physical  domain  is  chosen  as  one  quarter  of  an  ellipsoid-like  space. 

Details  about  the  coordinate  system,  governing  equations  and  boundary  conditions, 
and  the  numerical  method  are  provided  in  Reference  18  which  is  an  addendum  to  this 
report. 

The  drag,  lift,  and  moment  coefficients  are  evaluated  in  dimensional  form  as  follows. 


^  — n'n*e£)  dS'  -f  f  n^r'^eo  dS' 
s  Js 

(8) 

Fl  =  i 

(  —p'n-ei  dS'  -f-  f  n'T'-ec  dS' 
s  Js 

(9) 

M'  =  f  r' XT'  dS', 

(10) 

where  e^)  =  (— sinoji  -j-  cosaj),  ei  =  (cosati  -f  smaj).  S'  denotes  the  surface  of  the 
droplet,  n  is  the  outward  unit  normal  vector  at  the  surface,  r'  is  the  position  vector  from 
the  center  of  the  droplet,  and  t'  is  the  viscous  stress  tensor. 
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The  repelling  force  is  assumed  positive.  Due  to  the  geometrical  symmetry,  only  the 
z-component  of  moment  exists,  and  counter-clockwise  direction  is  assumed  positive. 

The  non-dimensional  coefficients  of  drag,  lift,  and  moment  are  defined  respectively  as 
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The  moving  droplets  change  their  directions  due  to  the  interactions  with  the  surround¬ 
ing  fluid.  However,  in  order  to  insure  its  accuracy,  the  computational  grid  must  satisfy 
the  condition  that  the  lift  and  torque  are  zero  for  a  single  droplet  moving  in  any  direction. 


3.2  Results  and  discussion 

We  first  tested  the  three-dimensional  code  by  solving  the  steady  axisymmetric  flow 
past  a  single  droplet  at  Reynolds  number  100.  Since  the  code  developed  solves  for  three 
Cartesian  components  of  velocity  in  the  transformed  grid,  an  axisymmetric  test  calculation 
still  exercises  the  fully  three-dimensional  aspects  of  the  code. 

Here,  we  discuss  the  flow  generated  by  an  impulsively  started  single  droplet  into  an 
initially  quiescent  fluid  using  the  three-dimensional  solution  procedure  for  viscosity  ratio 
of  25  and  a  density  ratio  of  300  (liquid  to  gas)  at  Reynolds  number  100.  The  time- 
dependent  solution  converges  asymptotically  to  the  steady-state  results.  Table  5  lists  the 
drag  coefficients  {Cdp  and  Cdv  are  respectively  the  pressure  and  viscous  parts  of  Cd) 
as  a  function  of  the  computational  grid  density  which  are  in  good  agreement  with  the 
correlation  of  Rivkind  &  Ryskin  [19].  Table  6  shows  two  angles  measured  from  the  front 
stagnation  point  where  is  the  angle  at  which  the  surface  vorticity  changes  its  sign  and 
62S  is  the  angle  at  which  the  surface  velocity  is  zero. 

The  calculations  were  performed  for  three  different  grids,  Ni  x  N2  x  N3  and  Nu  x 
N21  X  N^i  =  20  X  1 1  X  32  and  10  x  1 1  x  32,  30  x  15  x  48  and  15  x  15  x  48,  and  40  x  21  x  64 
and  20  X  21  X  64,  in  a  computational  domain  having  an  outer  boundary  located  at  21 
droplet  radii  from  the  droplet  center.  The  solutions  from  the  three  different  grids  were 
stable  and  smooth,  and  each  takes  dimensionless  times  of  20  to  reach  steady  state. 
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Tables  5  and  6  show  that  the  results  of  the  30  x  15  x  48  grid  differ  by  only  0.8  %  in 
the  drag  coefficient  and  0.7  %  in  the  separation  angle  from  those  of  the  40  x  21  x  64  grid. 
Thus,  for  computational  economy  (a  long  time  period  is  needed  for  our  time-dependent 
solution  of  the  moving  droplets),  we  selected  the  medium  size  grid  30  x  15  x  48  for  the 
external  flow  and  15  x  15  x  48  in  the  liquid  phase  for  the  remaining  calculations. 

For  the  interaction  of  two  moving  droplets,  calculations  were  performed  for  initial 
dimensionless  distance  dg  =  2  and  9  at  initial  Reynolds  number  100  for  the  time  period 
0  <  f  <  250.  This  period  is  about  one  quarter  of  the  lifetime  for  a  vaporizing  droplet 
injected  into  a  combustor,  and  so  major  changes  in  transport  processes  and  flow  field 
would  be  expected  to  occur.  The  viscosity  and  density  ratios  (liquid  to  gas)  for  the 
droplet  is  25  and  300,  respectively,  which  are  typiccd  of  a  liquid-hydrocarbon  fuel  in  a 
high-temperature,  high-pressure  surrounding  gas  generally  encountered  in  gas  turbine 
combustors. 

Figure  10(a)  shows  ihe  trajectory  of  one  droplet  for  =  2,  where  the  numbers  on 
the  line  denote  the  dimensionless  time  measured  from  injection.  The  final  location  of 
the  droplet  at  f  =  250  is  (dtjt)  =  (2.505,  212.6).  Now,  since  the  final  location  of  a 
single  droplet  at  f  =  250  in  the  case  of  no  drag  and  lift  forces  acting  on  it  would  be 
{dt,lt)  =  (2,  250),  the  figure  indicates  that  the  two  droplets  are  repelling  each  other  as 
well  as  decelerating,  and  the  change  of  distance  (normalized  by  droplet  radius)  is  much 
higher  (37.4  vs.  0.505)  due  to  the  drag  than  due  to  the  lift  force.  Figure  10(b)  shows  the 
trajectory  of  one  droplet  for  dg  =  9.  The  final  location  of  the  droplet  at  t  =  250  is  (dj.  It) 
=  (8.968,  213.1).  The  theoretical  final  location  of  a  single  droplet  at  <  =  250  in  the  case 
of  no  drag  and  lift  forces  would  be  (duh)  =  (9,  250);  therefore.  Figure  10(b)  indicates 
that  the  two  droplets  are  weakly  attracting  each  other  as  well  as  decelerating,  and  the 
change  of  distance  is  much  higher  due  to  the  drag  than  due  to  the  lift  force. 

Figure  11  shows  the  lift  coefficients  of  the  droplet  for  do  =  2  and  9  as  a  function  of 
the  instantaneous  Reynolds  number,  where  the  repelling  force  is  taken  as  positive.  The 
lift  coefficient  for  d,,  =  2  is  positive  during  the  time  period  0  <  t  <  250,  and  becomes 
gradually  smaller  in  time  because  the  distance  between  the  droplets  is  increasing  and  their 
directions  of  motion  are  changing  due  to  the  repelling  force.  On  the  other  hand,  the  lift 
coefficient  for  do  =  9  is  negative  during  that  time  period  but  slowly  goes  towards  zero. 
The  change  of  the  distance  between  the  two  droplets  for  do  =  9  case  is  negligibly  small 
as  shown  in  Figure  10(b).  Thus,  this  result  also  indicates  that  the  lift  coefficient  slowly 
approaches  zero  when  Reynolds  number  decreases  with  the  fixed  dimensionless  distance 
9.  The  figure  also  shows  that  rapid  change  occurs  initially  due  to  an  impulsive  start  of 
the  droplets  and  the  quasi-steady  state  occurs  when  t  >  30. 

Figure  12  shows  the  drag  coefficients  of  the  droplet  for  dg  —  2  and  9  as  a  function 
of  instantaneous  Reynolds  number.  The  drag  coefficient  for  dg  =  2  is  higher  than  that 
for  dg  =  9.  In  earlier  time,  the  difference  is  greater  but  becomes  gradually  smaller  as 
time  increases.  The  reason  is  that  the  droplets  for  dg  =  2  are  repelling  each  other  and 
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the  distance  between  them  is  increasing,  and  so  the  drag  becomes  gradually  smaller.  The 
drag  coeflBcient  for  do  =  9  is  almost  identical  to  that  for  a  single  droplet  (slightly  higher 
than  for  a  single  droplet  when  t  >  30).  We  found  in  Section  2  that  for  the  fixed  spheres 
and  do  >  4,  the  drag  coefficient  is  almost  identical  to  that  of  a  single  sphere.  Figure 
12  also  shows  large  drag  coefficient  initially  which  indicates  that  large  shear  stress  and 
pressure  occur  initially  due  to  impulsive  start. 

The  moment  coefficients  are  discussed  in  the  addendum.  That  paper  also  discusses  the 
sudden  drop  in  the  speed  at  the  moment  of  droplet  injection  and  the  monotonic  decrease 
with  time  following  injection.  The  speed  for  do  =  2  is  lower  than  do  =  9  due  to  the  larger 
drag  of  do  =  2. 


4.  Three-dimensional  flow  computation  for  two  interacting, 
vaporizing,  moving  droplets. 

We  consider  an  unsteady,  three-dimensional,  variable  property,  laminar  flow  generated 
by  two  identical  (spherical)  droplets  injected  into  an  initially  quiescent,  high-temperature, 
high-pressure  air  environment.  Then,  the  droplets  move  side  by  side  in  the  same  plane  as 
shown  in  Figure  9.  The  coordinate  system  is  exactly  same  as  that  used  in  Section  3. 

Since  our  goal  is  to  study  the  flow  interaction  with  the  two  droplets,  we  present  the 
equations  for  the  gas  and  liquid  phases.  The  governing  equations  in  both  phases  and 
the  boundary  conditions  are  nondimensionalized  using  the  initial  droplet  radius  aj,  as  the 
characteristic  length,  the  initial  droplet  velocity  U^  o  ^  the  characteristic  velocity,  and 
the  initial  air  temperature  ^  as  the  characteristic  temperature.  The  properties  in  the 
gas  phase  are  normalized  by  those  at  infinity,  and  the  properties  in  the  liquid  phase  are 
normalized  by  initial  droplet  properties. 


Gas  phase 


^  +  V.(;,V)  =  0 , 


dV,  dp\ 


2 


(U) 

(15) 


+  v.(pvr))  +  (h,  -  +  v-{pVY,)) 

^  V-(kVT)  +  -  K)pDVYi) , 


RePr 


ReSc 


(16) 
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(17) 


dpYj 


+  v*(/,vy})  = 


-V‘{pDVYj) , 


dt  ReSc 

where  V,  is  the  velocity  of  the  center  of  mass  of  the  two  droplets, 


is  the  deviatoric 


stress  tensor,  hj  and  ha  are  the  enthalpies  of  fuel  vapor  and  air,  respectively,  Yj  is  the 
fuel  mass  fraction,  and  D  is  the  diffusivity. 


Liquid  phase 


V-V,  =  0  , 

(18) 

(19) 

f)T  9 

(20) 

where  Vj  is  the  droplet  velocity  ajid  the  subscript  I  denotes  the  liquid  phase. 

The  governing  equations  are  transformed  to  the  generalized  coordinates  which 

allow  for  any  three-dimensional  body  of  arbitrary  shape.  The  numerical  integration 
of  the  equations  is  performed  using  a  computational  cubic  mesh  with  equal  spacing 
{8(  =  Srj  =  S(^  =  1).  The  conditions  at  the  interface  were  derived  from  principle  of  con¬ 
tinuity,  conservation,  and  thermodynamic  equilibrium.  The  overall  behavior  pf  droplet 
vaporization  is  strongly  dependent  on  the  evaluation  of  local  properties  at  the  interface 
where  complex  nonlinear  transport  phenomena  take  place.  A  quasi-linearization  method 
w<is  adopted  to  evaluate  iteratively  the  nonlinear  terms  in  the  interface  conditions.  Kim 
et  al.  [20]  will  provide  details  of  the  boundary  and  interface  conditions. 

The  code  we  have  developed  will  take  a  large  amount  of  cpu  hours  on  a  Cray  Y-MP 
due  to  the  integration  for  a  long  time  period  to  capture  an  unsteady  behavior  of  droplets 
and  three-dimensionality  of  the  problem.  To  check  our  numerical  algorithm,  we  solved  the 
flow  generated  by  a  single  cold  n-octane  fuel  droplet  suddenly  injected  into  a  hot  and  high 
pressure  gas  stream  by  using  a  coarse  grid  (20  x  11  x  32  in  the  gas  phase  and  12  x  1 1  x  32 
in  the  liquid  phase).  The  values  of  the  physical  parameters  used  are  given  in  Table  7. 

Figure  13  shows  the  history  of  the  drag  coefficients  for  a  single  vaporizing  droplet 
during  the  time  period  0  <  t  <  250.  It  is  shown  that  a  large  drag  coefficient  is  found 
initially  which  indicates  that  large  values  of  shear  stress  and  pressure  occur  initially  due  to 
the  impulsive  start.  During  the  initial  relaxation  period,  the  drag  coefficient  falls  rapidly. 
Subsequently,  the  drag  coefficient  decreases  slowly  due  to  the  blowing  effect  on  the  droplet 
surface  despite  a  reduction  in  the  Reynolds  number. 

Figure  14  shows  the  average  Nusselt  number  as  a  function  of  dimensionless  time.  The 
average  Nusselt  number  is  defined  as  follows. 
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(21) 


Nu„„  = 


2a[h' 


oo 


where  '  denotes  dimensional  quantity  and  k  is  the  thermal  conductivity  of  the  gas  phase 
normalized  by  the  thermal  conductivity  at  infinity.  Figure  14  shows  that  the  Nusselt 
number  is  decreasing  with  time.  The  Nusselt  number  is  the  dimensionless  rate  of  the  heat 
transfer,  and  the  decrease  of  the  Nusselt  number  is  mainly  attributed  to  the  reduction  of 
the  Reynolds  number  due  to  the  droplet  deceleration  under  the  action  of  the  drag. 

Figure  15  show  the  average  Sherwood  number  as  a  function  of  dimensionless  time. 
The  average  Sherwood  number  is  defined  in  dimensional  form  as  follows. 


= 


P'ooD'^ 


1 


2Tr(it(Yfoo 


j^pD[VYi)-jidS  , 


(22) 


where  '  denotes  dimensional  quantity  and  D  is  the  mass  diffusivity  of  the  gas  phase 
normalized  by  the  mass  diffusivity  at  infinity.  Figure  15  shows  that  the  Sherwood  number 
is  decreasing  with  time.  The  Sherwood  number  is  the  dimensionless  rate  of  the  mass 
transfer,  and  the  decrease  of  the  Sherwood  number  is  mainly  attributed  to  the  reduction 
of  the  Reynolds  number. 


5.  Axisymmetric  interactive  fuel  droplet  calculations 

This  axisymmetrical  calculation  considers  the  effects  of  variable  thermophysical  prop¬ 
erties,  transient  heating  and  internal  circulation  of  liquid,  deceleration  of  the  flow  due 
to  the  drag  of  the  droplet,  boundary-layer  blowing,  and  moving  interface  due  to  surface 
regression  as  well  as  relative  droplet  motion.  Details  of  the  analysis  are  provided  in  Ref¬ 
erence  [7]  which  is  included  as  an  addendum  to  this  report.  The  results  are  compared 
with  those  of  an  isolated  droplet  [21]  as  well  eis  those  of  the  two-droplet  system  [6]  to 
investigate  the  effect  of  the  presence  of  the  third  droplet.  The  interaction  effects  from  the 
downstream  or  upstream  droplet  are  identified.  Then  the  behavior  of  trailing  droplets, 
which  follow  the  first  two  or  three  droplets,  can  be  estimated  from  the  results  of  the 
first  two  or  three  droplets  on  account  of  the  nearly  periodical  nature  of  linear  droplet 
arrangements. 

The  variations  of  trajectory  with  time  and  drag  coeflficients  vs.  instantaneous  Reynolds 
number  for  the  cases  of  different  initial  spacings  are  illustrated  in  Figure  16.  D12  is  the 
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time-dependent  nondimensional  center-to-center  spacing  between  the  first  two  droplets 
and  D23  is  the  corresponding  spacing  between  the  second  and  third  droplets.  Results 
of  Ccise  1  with  large  initial  spacing  (D12=D23=12)  indicate  that  the  drag  coefficient  is 
smaller  for  the  second  droplet  than  for  the  lead  droplet  and  still  smaller  for  the  third 
droplet.  The  drag  coefficient  of  each  droplet  is  reduced  from  the  isolated  droplet  value. 
The  major  decrease  in  drag  occurs  in  the  first  two  droplets.  In  Case  2,  the  drag  coefficients 
of  the  lead  and  the  second  droplets  are  significantly  reduced  (comparing  curves  7  and  10, 
and  8  and  11,  respectively)  due  to  the  strong  interactions  when  the  first  two  droplets  are 
spaced  only  two  diameters  away.  Similar  trends  occur  for  the  downstream  droplet  pairs  of 
Case  3  (comparing  curves  8  and  14, and  9  and  15,  respectively).  The  third  droplet  in  Case 

2  has  a  higher  drag  coefficient  than  that  of  the  second  droplet  since  the  second  droplet 
is  better  shielded  by  the  droplet  before  it.  The  D23  thus  increases  with  time.  Also,  note 
that  the  drag  coefficient  of  the  third  droplet,  which  is  spaced  far  away  from  the  second 
droplet,  seems  to  be  approximately  independent  of  D12  as  indicated  in  curves  9  and  12. 
However,  it  strongly  depends  upon  D23  as  illustrated  in  curves  9,  12  and  15. 

Results  for  the  cases  of  sufficiently  large  spacing  (above  approximately  6  droplet  di¬ 
ameters)  show  that  the  flow  field  of  each  droplet  is  qualitatively  similar  to  that  of  an 
isolated  droplet,  although  the  transport  rates  are  reduced  along  the  downstream  direc¬ 
tion  due  to  the  cascade  effect  of  the  wake.  The  major  drops  in  transport  rates  occur  in 
the  first  two  droplets.  For  the  cases  of  small  droplet  spacings  (less  than  approximately 

3  droplet  diameters),  the  flow  field  of  downstream  droplets  can  be  significantly  altered 
due  to  the  interaction  effects  from  the  upstream  droplets.  Usually,  the  second  droplet  has 
the  lowest  drag  coefficient  since  it  interacts  with  both  neighboring  droplets.  However,  the 
difference  in  transport  rates  between  the  second  and  the  third  droplets  is  not  significant 
since  both  droplets  are  fully  protected  by  the  wake  of  the  first  droplet.  The  effect  of 
D23  on  the  behavior  of  the  first  droplet  is  insignificant.  However,  depending  upon  the 
values  of  D12  as  well  as  D23,  the  effect  of  D23  on  the  behavior  of  the  second  droplet  may 
become  significant.  The  general  qualitative  conclusions  drawn  from  the  two-droplet  study 
[6]  can  be  applied  to  two  neighboring  droplets  in  the  three-droplet  analysis.  The  correla¬ 
tions  for  heat  transfer  and  droplet  dynamics  have  been  developed  and  are  applicable  in  a 
one-dimensional  droplet  spray  calculation. 

In  addition  to  the  general  qualitative  conclusions  drawn  from  the  two-droplet  study, 
we  have  further  obtained  the  following  conclusions: 

(1)  Results  for  the  cases  of  sufficient!}'  large  spacing  (above  approximately  6  droplet 
diameters)  show  that  the  flow  field  of  each  droplet  is  qualitatively  similar  to  the  field  of  an 
isolated  droplet,  although  the  transport  rates  are  reduced  along  the  downstream  direction 
due  to  the  cascade  effect  of  the  wake. 

(2)  For  the  cases  of  small  droplet  spacings  (less  than  approximately  3  droplet  di¬ 
ameters),  the  flow  field  of  downstream  droplets  can  be  significantly  altered  due  to  the 
interaction  effects  from  the  upstream  droplets. 
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(3)  The  major  drop  in  transport  rates  occur  in  the  first  two  droplets.  Usually,  the 
recirculating  wake  from  the  lead  droplet  has  dissipated  a  significant  portion  of  the  mo¬ 
mentum  and  energy  upstream  of  the  second  droplet.  As  a  result,  the  further  reduction  in 
transport  rates  from  the  second  to  the  third  droplet  becomes  very  limited. 

(4)  Usually,  the  second  droplet  has  the  lowest  drag  coefficient  since  it  interacts  with 
both  neighboring  droplets.  However,  the  difference  in  transport  rates  between  the  second 
and  the  third  droplets  is  not  significant. 

(5)  The  effect  of  D23  on  the  behavior  of  the  first  droplet  is  insignificant.  However, 
depending  upon  the  values  of  D12  as  well  as  D23,  the  effect  of  D23  on  the  behavior  of  the 
second  droplet  is  not  negligible. 

(6)  The  correlations  for  heat  transfer  and  droplet  dynamics  have  been  developed  [7] 
and  are  applicable  to  droplet  spray  calculation. 


6.  LOX  vaporization  at  subcritical  conditions 

The  second  task  involves  the  axisymmetric  calculation  of  an  isolated  LOX  (liquid 
oxygen)  droplet.  A  detailed  analysis  of  LOX  droplet  heating  and  vaporization  and  the 
mixing  of  the  oxygen  and  fuel  vapors  in  a  high  temperature,  low  pressure  (10  -  20  atm), 
convective  environment  has  been  conducted.  This  problem  features  very  complex  property 
calculations  along  with  very  high  transfer  number  (  w  0(100) )  and  high  density  ratio.  The 
variation  of  liquid-phase  density  with  temperature  is  significant.  The  LOX  may  experience 
the  thermal  expansion  during  the  transient  heating  and  the  surface  moving  velocity  cannot 
be  neglected.  In  order  to  account  carefully  for  unsteadiness  due  to  density  variation,  the 
variable  liquid-density  with  primitive  variables  {Vr,Vz,pi  and  pi)  are  employed.  Note 
that  in  previous  codes  [21,6],  the  liquid-phase  momentum  equations  are  formulated  by  a 
stream  function-vorticity  approach.  The  modifications  of  the  code  to  adapt  the  primitive 
variables  for  the  liquid  phase  as  well  as  to  employ  more  accurate  interface  conditions  have 
been  developed. 

The  general  considerations  include  :  the  effects  of  variable  thermophysical  properties, 
real  gas  behavior,  transient  heating  and  internal  circulation  of  liquid,  deceleration  of  the 
flow  due  to  the  drag  of  the  droplet,  boundary-layer  blowing,  and  moving  interface.  An 
implicit  finite-difference  scheme  has  been  developed  to  solve  the  complete  set  of  Navier- 
Stokes,  energy,  and  species  equations.  Some  preliminary  results  have  been  presented  in 
Chiang  and  Sirignano  [22].  The  results  are  presented  for  the  low  pressure  case.  The 
interesting  phenomena  due  to  the  large  surface  blowing  and  surface  boiling  are  examined. 

We  have  learned  in  all  calculations  that  the  LOX  droplet  surface  very  quickly  reaches 
the  boiling  temperature,  which  is  about  120  K  for  the  case  of  10  atms  ambient  pressure. 
However,  the  droplet  core  remains  cold  since  the  conductivity  is  low  and,  at  the  early  time, 
the  circulation  strength  is  not  strong  enough  to  convect  much  energy.  All  the  available 
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heat  transfer  from  the  gas  phase  is  used  to  heat  the  surface  and  is  also  utilized  for  the 
latent  heat  of  vaporization.  The  surface  boiling  persists  throughout  the  droplet  lifetime. 
The  very  rich  and  cold  oxygen  vapor  surrounding  the  droplet  is  depicted  in  Figure  17. 
The  large  gradient  occurs  at  the  front  stagnant  region  as  we  expect.  The  LOX  droplet 
would  require  much  of  energy  for  the  heat  of  vaporization  rather  than  for  the  heating  the 
internal  fluid.  As  a  result,  the  transient  droplet  heating  is  responsible  for  the  unsteady 
thermal  behavior  of  the  LOX  droplet.  The  high  surface  blowing  velocity  has  a  significant 
effect  on  the  flow  structure  and  modifies  the  flow  separation  angle  and  wake  length,  etc. 
The  drag  coefficient,  Nusselt  and  Sherwood  numbers  are  reduced  to  values  below  their 
corresponding  stagnant  values.  A  copy  of  the  preprint  [22]  describing  this  research  is 
attached  as  an  addendum  to  this  report. 


7.  LOX  vaporization  at  near-critical  or  supercritical  conditions 

Most  rocket  engine  combustors  are  operated  in  the  very  high  pressure  domain  where 
LOX  droplets  may  vaporize  at  near  or  even  super-critical  conditions.  Hence,  the  cur¬ 
rent  objective  of  this  task  is  concentrated  on  an  understanding  of  fundamental  transport 
processes  underlying  high  pressure  LOX  droplet  vaporization. 

The  extension  of  the  present  low  pressure  model  to  cover  the  pressure  range  from  the 
subcritical  region  to  the  supercritical  region  is  currently  underway.  The  solubility  of  the 
fuel  vapor  in  the  liquid  phase  makes  a  multicomponent  droplet  formulation  necessary. 
Also,  a  new  technique  is  needed  to  compute  mass  fraction  and  thermophysical  properties 
at  gas/liquid  phase  equilibrium  and  to  account  for  the  droplet  regression  rate  due  to 
vaporization  as  well  as  phase  change  across  the  critical  surface,  with  continuous  density 
and  temperature  gradients  at  the  surface. 

Several  modifications  associated  with  physical  behavior  have  been  conducted  to  extend 
the  previous  low  pressure  LOX  droplet  vaporizing  model  to  deal  with  the  ambient  pressure 
in  the  range  from  the  near-critical  region  to  the  supercritical  region.  They  are  summarized 
below. 

1 .  A  comprehensive  calculation  of  vciriations  of  thermophysical  properties  with  respect 
to  temperature  for  each  species  component  is  performed.  The  high  pressure  correction, 
the  mixing  rules  to  evaluate  the  thermodynamic  properties,  and  the  quantum  gas  behavior 
of  hydrogen  vapor  are  considered.  The  gas-phase  and  liquid-phase  properties  tables  have 
been  constructed  in  order  to  evaluate  the  thermodynamic  properties  efficiently  by  a  two- 
variable  linear  interpolation. 

2.  The  real  gas  effect  in  the  high  pressure  case  is  taken  into  account  by  incorporating  a 
compressibility  factor  in  the  gas-phase  equation  of  state.  A  modified  Redlich-Kwong  equa¬ 
tion  of  state  with  the  mixing  rules  of  Chueh  and  Prausnitz  are  employed.  The  enthalpy 
of  vaporization  is  derived  from  the  equation  of  state  by  using  the  fugacity  equation. 


20 


Since  the  swelling  effect  of  the  LOX  vaporization  is  significant,  a  primitive  variable 
approach  is  formulated  to  analyze  the  internal  motion  and  thermal  and  mass  transport 
within  the  droplet. 

3.  The  solubility  of  the  fuel  vapor  in  the  liquid  phase  has  been  considered  by  a 
multicomponent  formulation.  The  principle  of  equal  temperature,  pressure,  and  fugacity 
of  both  phases  for  phase  equilibria  is  used  to  compute  the  surface  mass  fraction  of  species 
for  gas  and  liquid.  A  brief  description  of  the  new  phase  equilibrium  conditions  at  the 
droplet  interface  is  given  below. 

Thermodynamic  equilibrium  of  a  mixture  gives 

Ti  =  P*  =  P®,  fl  =  ff  (23) 

which  ultimately  presents  the  following  system  of  equations  to  be  solved: 

^02  X  P)  =  X  P) 

^  =  ^H2  ^  (^02 »  ^H2  » ^5  (24) 

(25) 

A  phase  equilibrium  diagram  and  enthalpy  of  vaporization  from  this  analysis  is  given 
in  Figure  18. 

The  single-component  interface  equations  are  replaced  by  the  following  interface  equa¬ 
tions  for  a  multicomponent  fluid. 


Conservation  of  species: 
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Conservation  of  energy: 
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(26) 


(27) 


where  m„  =  p{V  —  r)r+ 


4.  A  new  technique  is  being  tested  to  locate  the  gas/liquid  interface  and  to  compute 
the  regression  rate  due  to  vaporization  as  well  as  phase  change  across  the  critical  surface, 
with  continuous  density  and  temperature  gradients  at  the  surface. 
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A  calculation  is  performed  to  study  vaporization  behavior  before  the  droplet  reaches 
the  critical  mixing  state.  The  changes  of  properties  at  the  interface  occur  very  rapidly 
during  the  droplet  vaporization.  Diffusion  length  scales  are  reduced  at  the  elevated  pres¬ 
sure.  The  transport  rates  are  faster  than  those  in  the  low  pressure  case;  hence,  the  droplet 
lifetime  is  shorter.  As  a  result,  a  very  small  time  step  is  mandatory  to  ensure  accuracy. 
Figure  19  indicates  that  radius  increases  with  time  due  to  the  swelling  effect  of  droplet  as 
well  as  the  condensation  of  vapor.  The  reduced  enthalpy  of  vaporization  also  contributes 
to  the  fast  vaporization  of  the  droplet.  The  enthalpy  of  vaporization  shown  in  Figure  20 
is  significantly  reduced  from  the  conventional  “latent  heat” .  A  remarkable  condensation 
of  gas-phase  vapor  occurs  at  the  early  droplet  lifetime.  The  vaporization  begins  as  the 
surface  temperature  increases. 

Figure  21  presents  the  variations  of  compressibility  factors  (averaged  around  the  az¬ 
imuthal  direction)  for  the  gas  and  liquid  phases  at  the  interface.  When  the  two  com¬ 
pressibilities  are  equal,  the  critical  mixing  state  is  reached.  The  apparent  trend  of  two 
phases  merging  into  one  continuous  phase  is  shown.  Since  the  critical  mixing  state  highly 
depends  on  temperature  and  mass  fraction,  the  droplet  surface  at  different  locations  will 
reach  the  critical  mixing  state  at  different  times.  The  diffusivity  for  the  liquid  phase  is, 
in  general,  two  orders-of-magnitude  less  than  that  of  the  gas  phase;  as  a  result,  most  of 
the  diffusing  vapors  are  concentrated  in  the  boundary-layer  region  at  the  interface. 

Various  kinds  of  approaches  and  numerical  schemes  and  algorithms  are  under  exami¬ 
nation  to  address  the  following  difficulties  occurring  during  the  computation. 

1.  The  surface  temperature  rises  very  quickly  for  the  high  pressure  case.  Since  the 
liquid- vapor  equilibrium  conditions  require  a  separate  iterative-scheme  to  solve  for  surface 
mass  fraction  and  temperature,  we  have  temporarily  de-coupled  the  momentum  equations 
from  the  energy  and  species  conditions  and  solved  them  sequentially.  We  found  that  the 
present  method  tends  to  overestimate  the  blowing  velocity  at  the  early  droplet  time.  A 
relaxation  factor  is  hence  applied  to  decrease  artificially  the  vaporization  rate  during  the 
early  period. 

2.  An  appropriate  numerical  boundary  condition  at  the  droplet  center  must  be  devel¬ 
oped  for  the  present  cylindrical-coordinate  velocity  formulation  and  spherical-grid  system. 
Otherwise,  local  errors  occur. 

3.  Some  information  about  thermophysical  properties  at  the  critical  state  is  still  being 
sought. 

4.  The  entire  calculation  has  been  very  time-consuming.  Especially,  computation  of 
the  compressibility  factor  (which  requires  an  iterative  procedure  to  solve  a  cubic  polyno¬ 
mial)  utilizes  a  large  portion  of  computing  time.  We  are  pursuing  alternative  strategies 
(i.e.  by  a  linear  interpolation  method  to  fetch  the  values  from  existing  tables;  however,  it 
would  require  huge  memory  storages)  for  simplification  of  this  part  of  the  calculation. 

Currently,  we  have  advanced  our  high-pressure  LOX  model  to  deal  with  the  vapor¬ 
ization  at  supercritical  conditions.  The  major  difficulties  in  the  calculation  are  :  1)  the 
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formulation  of  boundary  conditions,  especially  to  trace  the  location  of  the  “droplet  inter¬ 
face”,  2)  the  property  calculation  of  the  “gas/liquid  puff”  at  the  interface.  Basically,  we 
have  developed  two  modelling  approaches,  which  are  currently  undergoing  investigation, 
to  handle  the  interface  conditions.  Here,  we  briefly  describe  our  strategies  step-by-step 
to  make  calculations  at  supercritical  conditions. 

(1)  Determine  the  critical  mixing  temperature  with  respect  to  various  pressure  levels. 
The  Redlich-Kwong  equation  of  state  with  the  mixing  rules  of  Chueh  and  Prausnitz  is 
employed  for  vapor-liquid  equilibria  calculation  for  H2IO2  system.  The  critical  mixing 
temperature  is  determined  when  the  compressibilities  (or  mass  fractions  of  a  component) 
at  both  phases  reach  the  same  value.  Figure  22  demonstrates  that  the  critical  mixing  tem¬ 
perature  and  oxygen  mass  fraction  at  the  interface  reduce,  while  compressibility  increcises 
as  pressure  increases. 

(2)  The  droplet  “pseudo- interface  location”  is  determined  by  the  critical  mixing  tem¬ 
perature.  Note  that  the  pseudo- interface  will  slightly  deviate  from  its  original  spherical 
shape  due  to  the  convective  effects  in  both  gas  and  liquid  phases.  In  calculating  the 
vaporization  rate,  it  is  still  assumed  that  the  droplet  remains  spherical  (by  the  use  of  an 
averaged  radius). 

(3)  The  first  approach  assumes  that  there  is  no  gas/liquid  interface.  The  flowfield 
is  continuous  from  the  gas  phase  to  the  liquid  phase  and  the  problem  is  solved  as  a 
single-phase  problem.  The  solutions  strongly  depend  on  the  property  evaluation  at  the 
interface. 

(4)  Since  there  is  limited  availability  of  the  critical  point  thermo-physical  properties 
data  in  the  literature,  we  assume  that  all  thermo-physical  properties  are  also  continuous 
from  the  liquid  to  the  gas  phase.  The  consideration  of  singular  behavior  of  some  thermo¬ 
physical  properties  is  excluded  in  the  present  study.  A  second-order  linear  interpolation 
scheme  is  employed  to  obtain  the  values  at  critical  mixing  state.  However,  we  have 
learned  that  the  direct  calculation  of  density  using  this  approach  will  lead  to  a  failure  in 
the  calculation  of  vaporization  rate.  Hence,  the  density  is  obtained  from  compressibility 
relations.  The  variation  of  mass  with  respect  to  time  is  shown  in  Figure  23.  After  an 
early  condensing  period,  the  droplet  reaches  the  critical  mixing  state  and  the  evaporation 
rate  increases  dramatically  beyond  the  critical  point. 

(5)  Another  modelling  approach  assumes  that  there  exists  a  spherical  “pseudo-interface.” 
The  gas  and  liquid  phase  equations  are  solved  separately,  and  the  solutions  are  matched 
at  the  interface.  In  this  approach,  constant-temperature  (at  its  critical  mixing  state)  and 
constant- mass-fraction  interface  conditions  are  employed.  The  shear  (or  normal)  stress 
conditions  are  applied  even  though  the  surface  tension  vanishes  at  the  interface.  We  still 
suffered  some  numerical  difficulties  in  this  approach  and  the  preliminary  numerical  results 
are  not  comparable  to  the  existing  data. 

In  parallel  to  the  high-pressure  LOX  calculation,  we  are  modifying  the  previous  pro¬ 
gram  for  the  interactive  hydrocarbon-fuel  droplet  to  take  LOX  or  hydrogen  droplet  va- 
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porization  into  account.  The  modifications  include:  (1)  a  primitive  variables  approach  for 
the  liquid  phase,  (2)  thermodynamics  property  evaluations  for  oxygen  and  hydrogen  fuel, 
and  (3)  more  accurate  formulation  of  the  interface  to  include  droplet  surface  expansion. 

We  currently  are  addressing  several  technical  difficulties  listed  below: 

1.  There  is  a  need  for  a  suitable  initial  variable  profile  to  start  the  calculation.  The 
pressure  field  (and  mass  residual  sum)  between  two  droplets  could  destroy  the  solutions 
since  a  large  magnitude  of  vaporization  occurs  if  droplets  suddenly  encounter  an  inappro¬ 
priate  temperature  profile  at  the  beginning  of  the  simulation. 

2.  The  computation  will  be  very  costly.  The  whole  primitive  variables  calculation 
for  a  vaporizing  LOX  or  hydrogen  droplet  requires  a  very  small  time  step  to  guarantee 
a  converged  solution  at  the  droplet  interface.  In  addition,  the  grid  generation  routine 
requires  more  iterations  to  keep  track  of  the  rapidly  recessing  droplet  surface. 

We  will  modify  the  current  computational  model  to  improve  the  numerical  stability 
problem.  Also,  a  simple  formulation  of  droplet  vaporization  may  be  constructed  in  order 
to  conserve  our  very  limited  CRAY  Y-MP  — .  More  detailed  thermophysical  property 
information  for  a  hydrogen  droplet  is  ci  -rcmly  being  compiled. 


Acknowl'^dgement 

We  would  like  to  thank  Mr.  Lyle  Wiedeman  and  Dr.  Allen  Schiano  in  the  Office  of 
Academic  Computing  of  UCI  for  their  assistance  in  using  a  three-dimensional  graphic 
package  Application  Visualization  System  (AVS).  The  support  of  the  Pittsburgh  Super¬ 
computing  Center  and  the  San  Diego  Supercomputer  Center  under  a  block  grant  of 
Office  of  Academic  Computing  of  UCI  are  gratefully  acknowledged. 


24 


References 


1.  Sirignano,  W.  A.,  1983,  Fuel  droplet  vaporization  and  spray  combustion.  Prog. 
Energy  and  Comb.  Set.  291-322. 

2.  Patnaik,  G.,  1986,  A  numerical  solution  of  droplet  vaporization  with  convection. 
Ph.D.  Dissertation,  Carnegie-Mellon  University. 

3.  Tal(thau),  R.,  Lee,  D.  N.  &  Sirignano,  W.  A.,  1983,  Hydrodynamics  and  heat 
transfer  in  sphere  ci.3semblages  -  cylindrical  cell  models.  Int.  J.  Heat  Mass  Transfer 
26,  No.  9,  1265-1273. 

4.  Tal(thau),  R.,  Lee,  D.  N.  &  Sirignano,  W.  A.,  1984,  Heat  and  momentum  transfer 
around  a  pair  of  spheres  in  viscous  flow.  Int.  J.  Heat  Mass  Transfer  27,  No.  11, 
1253-1262. 

5.  Raju,  M.  S.  Sz  Sirignano,  W.  A.,  1990,  Interaction  between  two  vaporizing  droplets 
in  an  intermediate- Reynolds-number  flow.  Phys.  Fluids  A, 2(10),  1780-1796. 

6.  Chiang,  C.  H.  &  Sirignano,  VV.  A.,  1993a,  Numerical  analysis  of  interacting,  convect- 
ing,  vaporizing  fuel  droplets  with  variable  properties.  Int.  J.  Heat  Mass  Transfer, 
in  press, 

7.  Chiang,  C.  H.  ic  Sirignano,  W.  A.,  1993b,  Axisymmetric  calculations  of  thfee-droplet 
interactions.  Atomization  and  Sprays,  in  press. 

8.  Anderson,  D.  A.,  Tannehill,  J.  C.  &  Pletcher,  R.  H.,  1984,  Computational  Fluid 
Mechanics  and  Heat  Transfer.  Hemisphere  Publishing. 

9.  Kim,  L,  Elghobashi  S.  &  Sirignano,  W.  A.,  1993,  Three-dimensional  flow  over  two 
spheres  placed  side  by  side.  J.  Fluid  Mech.  246,  465-488. 

10.  Roos,  F.  W.  &:  Willmarth,  W,  W.,  1971,  Some  experimental  results  on  sphere  and 
disk  drag.  AIAA  J.  9,  285-291. 

11.  Clift,  R.,  Grace,  J.  R.  &  Weber,  M.  E.,  1978,  Bubbles,  Drops,  and  Particles.  .Aca¬ 
demic  Press,  New  York. 

12.  Taneda,  S.,  1956,  Experimental  investigation  of  the  wake  behind  a  sphere  at  low 
Reynolds  number.  J.  Phys.  Soc.  Japan  11  1104-1108. 

13.  Rosfjord,  T.  J.,  1974,  Experimental  and  theoretical  investigations  of  the  recirculat¬ 
ing  flow  region  between  two-dimensional,  parallel,  separated  jets.  Ph.D.  Disserta¬ 
tion,  Princeton  University. 


25 


14.  Dandy,  D.  S.  &  Dwyer,  H.  A.,  1990,  A  sphere  in  shear  flow  a  finite  Reynolds  number; 
effect  of  shear  on  particle  lift,  drag,  and  heat  transfer.  J.  Fluid  Mech.  216,  381-410. 

15.  Goldburg,  A.  &  Florsheim,  B.  H.,  1966,  Transition  and  Strouhal  number  for  the 
incompressible  wake  of  various  bodies.  Phys.  Fluids  9,  45-50. 

16.  Nakamura  I.,  1976,  Steady  wake  behind  a  sphere.  Phys.  Fluids  19,  5-8. 

17.  Kim,  I.  Pearlstein,  A.  J.,  1990,  Stability  of  the  flow  past  a  sphere.  J.  Fluid  Mech. 
211,  73-93. 

18.  Kim,  I.,  Elghobashi  S.  &  Sirignano,  W.  A.,  1992,  Three-dimensional  flow  computa¬ 
tion  for  two  interacting,  moving  droplets.  AIAA  preprint  92-0343. 

19.  Rivkind  V.  Y.  Sz  Ryskin,  G.,  1976,  Flow  structure  in  motion  of  a  spherical  drop  in 
a  fluid  medium  at  intermediate  Reynolds  numbers.  Fluid  Dyn.  11,  5-12. 

20.  Kim,  I.,  Elghobashi  S.  &  Sirignano,  W.  A.,  1993,  Three-dimensional  flow  computa¬ 
tion  for  two  interacting,  vaporizing,  moving  droplets  with  variable  properties,  to  be 
presented  at  AIAA  24th  Fluid  Dynamics  Conference,  Orlando,  Florida. 

21.  Chiang,  C.  H.,  Raju,  M.  S.  Sirignano,  W.  A.,  1992,  Numerical  analysis  ot  convect- 
ing,  vaporizing  fuel  droplet  with  variable  properties.  Int.  J.  Heat  Mass  Transfer, 
35,  No.  5,  pp.  1307-1324. 

22.  Chiang,  C.  H.  &  Sirignano,  W.  A.,  1991,  Axisymmetric  vaporizing  oxygen  droplet 
computations.  AIAA  preprint  9T0281,  presented  at  the  AIAA  29th  Aerospace  Sci¬ 
ences  Meeting,  Reno,  Nevada. 


26 


Professional  Personnel 


W.  A.  Sirignajio,  Professor,  Principal  investigator. 

S.  E.  Elghobashi,  Professor. 

I.  Kim,  Postdoctoral  research  associate,  Ph.  D.  January  1990  -  present. 

C.  H.  Chiang,  Postdoctoral  research  associate,  Ph.  D.  September  1990  -  present, 

Ph.  D.  candidate  prior  to  September  1990. 

Dissertation  Title:  Isolated  and  interacting,  vaporizing  fuel  droplets:  field  calculation 
with  variable  properties. 


Publications 

The  following  papers  and  one  book  article  have  resulted  at  legist  partially  from  the 
research  performed  under  this  research  program. 

1.  Raju,  M.  S.  &  Sirignano,  W.  A.,  1990,  Interaction  between  two  vaporizing  droplets 
in  an  intermediate- Reynolds-number  flow.  Phys.  Fluids  A, 2(10),  1780-1796. 

2.  Chiang,  C.  H.  &:  Sirignano,  W.  A.,  1991,  Axisymmetric  vaporizing  oxygen  droplet 
computations.  AIAA  preprint  91-0281,  presented  at  the  AIAA  29th  Aerospace  Sci¬ 
ences  Meeting,  Reno,  Nevada. 

3.  Sirignano,  W.  A.,  Chiang,  C.  H.,  Kim,  I.  k.  Elghobashi,  S.  E.,  1991,  Aerodynamic 
interactions  amongst  neighboring  droplets,  presented  at  the  4th  international  Sym¬ 
posium  on  Computational  Fluid  Dynamics. 

4.  Sirignano,  W.  A.,  1991,  Computational  challenges  in  spray  combustion,  presented 
at  the  U.  S.  -  Japan  Heat  Transfer  Seminar. 

5.  Delplanque,  J.  P.,  Chiang,  C.  H.  k  Sirignano,  W.  A.,  1991,  Numerical  simulation 
and  modelling  of  LOX  droplet  vaporization  at  supercritical  conditions,  presented  at 
the  28th  JANNAF  Combustion  Meeting. 

6.  Kim,  I.,  Elghobashi  S.  k  Sirignano,  W.  A.,  1991,  Three-dimensional  flow  interac¬ 
tions  between  two  neighboring  spheres,  presented  at  the  44^^  Annual  Meeting  of  the 
Division  of  Fluid  Dynamics  of  the  American  Physical  Society,  Scottsdale,  Arizona. 

7.  Kim,  I.,  Elghobashi  S.  k  Sirignano,  W.  A.,  1992,  Three-dimensional  flow  computa¬ 
tion  for  two  interacting,  moving  droplets.  AIAA  preprint  92-0343,  presented  at  the 
AIAA  30th  Aerospace  Sciences  Meeting,  Reno,  Nevada. 


27 


8.  Chiang,  C.  H.,  Raju,  M.  S.  &  Sirignano,  W.  A.,  1992,  Numerical  analysis  of  convect- 
ing,  vaporizing  fuel  droplet  with  variable  properties.  Int.  J.  Heat  Mass  Transfer, 
35,  No.  5,  pp.  1307-1324. 

9.  Sirignano,  W.  A.,  1993,  Computational  spray  combustion,  in  Numerical  Modelling 
in  Combustion,  (Chung,  T.  J.  ed.),  Hemisphere  Publication,  in  press. 

10.  Chiang,  C.  H.  &  Sirignano,  W.  A.,  1993a,  Numerical  analysis  of  interacting,  convect- 
ing,  vaporizing  fuel  droplets  with  variable  properties.  Int.  J.  Heat  Mass  Transfer, 
in  press.  Also  see  Preprint  No.  90-0357,  AIAA  28th  Aerospace  Sciences  Meeting. 

11.  Chiang,  C.  H.  &  Sirignano,  W.  A.,  1993b,  Axisymmetric  calculations  of  three- 
droplet  interactions.  Atomization  and  Sprays,  in  press.  Also  presented  at  the  Fifth 
international  ICLASS  Conference  on  Liquid  Atomization  and  Spray  Systems. 

12.  Kim,  I.,  Elghobashi  S.  k  Sirignano,  W.  A.,  1993,  Three-dimensional  flow  over  two 
spheres  placed  side  by  side.  J.  Fluid  Mech.  246,  465-488. 

13.  Kim,  I.,  Elghobashi  S.  &:  Sirignano,  W.  A.,  1993,  Three-dimensional  flow  computa¬ 
tion  for  two  interacting,  vaporizing,  moving  droplets  with  variable  properties,  to  be 
presented  at  the  AIAA  24th  Fluid  Dynamics  Conference,  Orlando,  Florida. 

14.  Sirignano,  W.  A.,  1993,  Fluid  dynamics  of  sprays  (invited  Freeman  Scholar  paper). 
J.  Fluids  Engineering,  in  press. 

15.  Sirignano,  W.  A.,  Delplanque,  J.  P.  h  Chiang,  C.  H.,  1993,  Liquid  propellant  droplet 
vaporization,  Driving  mechanism  for  rocket  combustion  instability,  to  be  presented 
at  the  First  International  Symposium  on  Liquid  Rocket  Combustion  Instability. 


Presentations 

1.  Chiang,  C.  H.,  “Numerical  analysis  of  interacting,  convecting,  vaporizing  fuel  droplets 
with  variable  properties,”  presented  at  the  AIAA  28th  Aerospace  Sciences  Meeting, 
Reno,  Nevada,  Jan.  1990. 

2.  Sirignano,  W.  A.,  “Combustion  instability  computational  analysis,”  presented  at  the 
JANNAF  Workshop  on  Numerical  Methods  in  Combustion  Instability,  Orlando,  FL, 
Feb.  1990. 

3.  Sirignano,  W.  A.,  “Fundamental  studies  of  droplet  interactions  in  dense  sprays,” 
presented  at  the  AFOSR  Contractors’  Meeting,  Atlanta,  GA,  June  1990. 


28 


4.  Sirignano,  W.  A.,  “Heat  and  mass  transfer  in  liquid-gas  spray  systems,”  presented 
at  the  ASME  Winter  Annual  Meeting,  Dallas,  TX,  Nov.  1990. 

5.  Chiang,  C.  H.,  “Axisymmetric  vaporizing  oxygen  droplet  computations,”  presented 
at  the  AIAA  29th  Aerospace  Sciences  Meeting,  Reno,  Nevada,  Jan.  1991. 

6.  Sirignano,  W.  A.,  ”Interacting  LOX  droplets  in  dense  sprays,”  presented  at  the 
AFOSR  Contractors’  Meeting,  Boulder,  CO,  June  1991. 

7.  Sirignano,  W.  A.,  “Aerodynamic  interactions  amongst  neighboring  droplets,”  pre¬ 
sented  at  the  4th  international  Symposium  on  Computational  Fluid  Dynamics,  UC 
Davis,  CA,  Sep.  1991. 

8.  Sirignano,  W.  A.,  “Computational  challenges  in  spray  combustion,”  presented  at 
the  U.  S.  -  Japan  Heat  Transfer  Seminar,  Oiso,  Japan,  Oct.  1991. 

9.  Sirignano,  W.  A.,  “Numerical  simulation  and  modelling  of  LOX  droplet  vaporization 
at  supercritical  conditions,”  presented  at  the  28th  JANNAF  Combustion  Meeting, 
Brooks  Air  Force  Base,  San  Antonio,  TX,  Oct.  1991. 

10.  Sirignano,  W.  A.,  “Axisymmetric  calculations  of  three-droplet  interactions,”  pre¬ 
sented  at  the  Fifth  international  ICLASS  Conference  on  Liquid  Atomization  and 
Spray  Systems,  Gaithersburg,  MD,  July  1991. 

11.  Kim,  L,  "Three-dimensional  flow  interactions  between  two  neighboring  spheres,” 
presented  at  the  44th  Annual  Meeting  of  the  Division  of  Fluid  Dynamics  of  the 
American  Physical  Society,  Scottsdale,  Arizona,  Nov.  1991. 

12.  Kim,  I.,  "Three-dimensional  flow  computation  for  two  interacting,  moving  droplets,” 
presented  at  the  AIAA  30th  Aerospace  Sciences  Meeting,  Reno,  Nevada,  Jan.  1992. 

13.  Sirignano,  W.  A.,  “Recent  Advances  in  fuel  droplet  and  spray  theory,”  UC  San 
Diego,  CA,  Jan.  1992. 

14.  Sirignano,  W.  A.,  “Recent  Advances  in  fuel  droplet  and  spray  theory,”  University 
of  Southern  California,  CA,  March  1992. 

15.  Sirignano,  W.  A.,  “Fundamentals  of  spray  combustion,”  AFOSR  Propulsion  Con¬ 
tractors  Meeting,  LaJolla,  CA,  June  1992. 

16.  Sirignano,  W.  A.,  “Fluid  dynamics  of  sprays  (invited  Freeman  Scholar  paper),” 
presented  at  the  ASME  Winter  Annual  Meeting,  Anaheim,  CA,  Nov.  1992. 


29 


17.  Kim,  I.,  “Three-dimensional  flow  computation  for  two  interacting,  vaporizing,  mov¬ 
ing  droplets  with  variable  properties,”  to  be  presented  at  the  AIAA  24th  Fluid 
Dynamics  Conference,  Orlando,  Florida,  July  1993. 


Interactions 

Many  of  the  papers  hsted  under  publications  were  presented  as  conference  papers. 
Even  journal  papers  usually  were  presented  in  a  precursor  form.  We  have  participated  in 
AFOSR  Contractors  Meetings  and  Workshops  as  requested  by  the  sponsor.  W.  A.  Sirig- 
nano  has  served  in  several  advisory  or  consultant  roles  with  NASA,  the  National  Research 
Council,  and  major  government  contractors.  In  particular,  he  has  served  as:  Member, 
NASA  Space  Science  and  Application  Advisory  Committee;  Member,  National  Research 
Council  Space  Studies  Board  and  Chairman,  Committee  on  Microgravity  Research;  Con¬ 
sultant,  United  Technologies  Research  Center;  Consultant,  Alhson  Gas  Turbine  Division 
of  General  Motors;  Consultant,  Sandia  National  Laboratories;  Member,  USRA  Science 
Council  for  Microgravity  Science  and  Applications.  Dr.  Sirignano  has  been  involved  in 
other  research  programs  funded  by  AFWAL,  AFAL,  and  ONR  that  have  some  relation¬ 
ship  to  the  research  performed  in  this  program.  Various  invited  seminars,  invited  papers 
or  articles,  and  honors  have  resulted  for  the  principal  investigator  based,  at  least  par¬ 
tially,  on  this  AFOSR  research.  Examples  are  the  1992  ASME  Freeman  Scholar  Award, 
which  involves  a  special  ASME  lecture  and  journal  article,  the  1992  AIAA  Combustion 
and  Propellants  Award,  upgrade  to  Fellow  of  the  AAAS,  and  the  1991  AIAA  Pendray 
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Nx  X  jVj  X  A3 

Cdp 

Cov 

Cd 

20  X  21  X  21 

0.706 

0.951 

1.657 

30  X  31  X  31 

0.683 

0.9.34 

1.617 

40  X  41  X  41 

0.676 

0.929 

1.605 

1.58 

Table  1.  Drag  coefficients  as  a  function  of  grid  density  at  Re  =  50, 
where  *  denotes  the  data  from  Roos  &:  Willmarth  (1971) 
and  also  Clift  et  ai  (1978). 


iVi  X  iVj  X  N3  Pof  Pot  K 

20  X  21  X  21  0.611  -0.1088  136.53 

30  X  31  X  31  0.606  -0.0982  138.10 

40  X  41  X  41  0.604  -0.0954  138.63  139.3 

Table  2.  Pressure  at  the  front  and  rear  stagnation  points  and  the 
separation  angle  measured  from  the  front  stagnation  point 
as  a  function  of  grid  density  at  Re  =  50,  where  *  denotes 
the  data  from  Clift  et  al.  (1978). 


Ni  X  jV2  X  .V3 

Cdp 

Cdv 

Cd 

20  X  21  X  21 

0.555 

0.593 

1.148 

30  X  31  X  31 

0.532 

0.582 

1.114 

40  X  41  X  41 

0.524 

0.581 

1.105 

1.09 

Table  3.  Drag  coefficients  as  a  function  of  grid  density  at  Re  =100, 
where  *  denotes  the  data  from  Roos  &  Willmaxth  (1971) 
and  also  Clift  et  al.  (1978). 


Nx  X  N2  X  Nz  Ppi  Pot  Os 

20  X  21  X  21  0.555  -0.0924  124.24 

30  X  31  X  31  0.555  -0.0819  125.74 

40  X  41  X  41  0.554  -0.0789  126.25  126.5 

Table  4.  Pressure  at  the  front  and  rear  stagnation  points  and  the 
separation  angle  measured  from  the  front  stagnation  point 
as  a  function  of  grid  density  at  Re  =  100,  where  *  denotes 
the  data  from  Taneda  (1956)  and  Clift  et  al.  (1978). 


Gas 


Liquid 


iVi  X  N2  X  N3 

X 

X 

Cdp 

Cdv 

Cd 

20  X  11  X  32 

10  X  11  X  32 

0.542 

0.586 

1.128 

30  X  15  X  48 

15  X  15  X  48 

0.520 

0.573 

1.092 

40  X  21  X  64 

20  X  21  X  64 

0.511 

0.571 

1.081 

1.08 

Table  5.  Drag  coefficients  as  a  function  of  grid  density  at  Re  =  100, 
where  *  denotes  the  data  from  the  correlation, 
of  Rivkind  &  Ryskin  (1976) 


Gas 

Liquid 

N1XN2X  N3 

;Vi,  X  N21  X  N31 

9u 

02, 

20  X  11  X  32 

10  X  11  X  32 

126.81 

146.50 

30  X  15  X  48 

15  X  15  X  48 

127.70 

150.13 

40  X  21  X  64 

20  X  21  X  64 

127.74 

151.20 

Table  6.  Angles  at  which  the  surface  vorticity  changes  its  sign  and 
the  surface  velocity  is  zero  as  a  function  of  grid  density 
at  Re  =  100,  where  the  angles  are  measured  from 
the  front  stagnation  point. 


Parameter  Value 

Initial  Reynolds  number  in  gas  phase  100 

Relative  velocity  of  droplet,  m/s  25 

Free  stream  temperature,  K  1250 

Combustor  pressure,  atm  10 

Prandtl  number  in  gas  phase  0.74 

Prandtl  number  in  liquid  phase  8.53 

Schmidt  number  in  gas  phase  2.33 

Molecular  weight  of  oxidizer,  kg/kmol  29.0 

1'  ’^  ]'  .ular  weight  of  fuel  (n-octane),  kg/kmol  114.2 
Vroplet  initial  temperature,  K  300 

Viscosity  ratio  (liquid  to  gas)  10.5 

Density  ratio  (liquid  to  gas)  247.9 

Specific  heat  ratio  (liquid  to  gas)  1.89 

Latent  heat,  cal/g  72.5 


Table  7.  Physical  parameters  used  in  the  calculation 


Figure  2.  Sketch  of  typical  streamlines  over  one  of  the  two  solid  spheres 

in  the  principal  plane  (x-z  plane)  at  Re  =  100  for  (a)  axisymmetric  flow 
(b)  dg  —  2;  (c)  dg  =  1.5. 


Figure  3(b)  Viscous  contribution  to  the  lift  coefficient  of  the  soUd  sphere 
as  a  function  of  at  i?e  =  50,  100,  and  150 


Figure  3(c)  Pressure  contribution  to  the  lift  coefficient  of  the  solid  sphere 
as  a  function  of  do  at  Re  =  50,  100,  and  150 
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Figure  6.  Velocity  vector  fields  in  the  principal  plane  for  a  liquid  sphere 
at  Re  —  100  for  do  =  1.5  :  (a)  external  flow;  (b)  internal  flow 


Figure  7.  Total  lift  coefficients  of  the  liquid  sphere  as  a  function  of 
do  at  Re  =  50,  100,  and  150 


Figure  12.  Drag  coefficients  as  a  function  of  instantaneous  Reynolds  number 
for  do  =  2  and  9,  where  the  numbers  on  the  curves  denote 
the  dimensionless  time. 


C-C  SPACING  BETWEEN  TWO  DROPLETS 


Legends: 

Case  1:  Initial  Ri  =  i,i?2  =  l,i?3  =  1,  D12=12,  D23=12,  Re=100 
Case  2:  Initial  =l,R2  =  l,i?3  =  1,  D12=4,  D23=12,  Re=100 
Case  3:  Initial  R^  =  l,i?2  =  l,i?3  =  1,  D12=12,  D23=4,  Re=100 

I.  D12,  Case  1;  2.  D23,  Case  1;  3.  D12,  Case  2; 

4.  D23,  Case  2;  5.  D12,  Case  3;  6.  D23,  Case  3; 

7.  Cd  for  the  Lead  Droplet,  Case  1;  8.  Cd  for  the  Second  Droplet,  Case  1; 

9.  Cd  for  the  Third  Droplet,  Case  1;  10.  Cd  for  the  Lead  Droplet,  Case  2; 

II.  Cd  for  the  Second  Droplet,  Case  2;  12.  Cd  for  the  Third  Droplet,  Case  2; 
13.  Cd  for  the  Lead  Droplet,  Case  3;  14.  Cd  for  the  Second  Droplet,  Case  3; 
15.  Cd  for  the  Third  Droplet,  Case  3;  16.  Cd  for  an  Isolated  Droplet; 
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Figure  16.  Time  variation  of  droplet  drag  coefficients  and  droplet  spacings 
for  the  cases  of  large  initial  droplet  spacings. 
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Time  =  5.00  ,  Instantaneous  Reynolds  Number  =  83.85 


Figure  17.  O-xygen  mass-fraction  contours  and  isotherms. 
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Figure  18.  Hydrogen/Oxygen  vapor-liquid  equilibria  and  enthalpy  of  vaporization 
for  liquid  oxygen  at  different  reduced  pressures. 


Change 


LOX  in  H2 


Figure  19.  Time  variations  of  droplet  radius,  m^lss  and  liquid  density. 


Enthalpy  of  Vaporization 


Enthalpy  of  Vaporization 
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Gas-Phase  Diffusion  Time 

Figure  20.  Time  variations  of  enthalphy  of  vaporization  for  two  components. 


Compressibilities  for  Gas  and  Liquid 


Gas-Phase  Diffusion  Time 


Figure  21.  Time  variation  of  compressibility  of  gas  and  liquid  phases. 
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Figure  22.  Variations  of  critical  mixing  temperature,  oxygen  meiss  fraction, 
and  compressibility  with  respect  to  pressure. 


Figure  23.  variations  of  droplet  radius  with  respect  to  time  for  the  cases  of 
two  different  pressure  levels. 
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Three-dimensional  flow  over  two  spheres  placed 

side  by  side 
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(Received  23  August  1991  and  in  revised  form  2S  July  1992) 

Three-dimensional  flow  over  two  identical  (solid  or  liquid)  spheres  which  are  held 
filed  relative  to  each  other  with  the  line  connecting  their  centres  normal  to  a  uniform 
stream  is  investigated  numerically  at  Reynolds  numbers  50,  100,  and  150.  We 
consider  the  lift,  moment,  and  drag  coefiBcients  on  the  spheres  and  investigate  their 
dependence  on  the  distance  between  the  two  spheres.  The  computations  show  that, 
for  a  given  Reynolds  number,  the  two  spheres  are  repelled  when  the  spacing  is  of  the 
order  of  the  diameter  but  are  weakly  attracted  at  intermediate  separation  distances. 
For  small  spacing,  the  vortical  structure  of  the  near  wake  is  significantly  different 
from  that  of  the  axisymmetric  wake  that  estabUshes  at  large  separations.  The 
partially  confined  flow  passing  between  the  two  spheres  entrains  the  flows  coming 
around  their  other  sides.  Our  results  agree  with  available  experimental  and 
numerical  data. 


1.  Introduction 

Flow  past  droplets  and  solid  particles  are  important  in  many  natural  and 
engineering  applications  such  as  air  pollution,  combustion  systems,  and  chemical 
processes.  Many  investigations  have  considered  the  interactions  between  droplets  or 
particles  and  the  simrounding  fluid  by  analytical  and  numerical  methods.  For' 
sufficiently  low  or  high  Reynolds  numbers,  a  theoretical  analysis  can  be  performed 
using  singular  perturbation  expansions  which  involve  linearization  or  the  boundary- 
layer  approximation.  For  flows  at  intermediate  Reynolds  numbers,  which  are  most 
common  in  engineering  applications,  it  is  necessary  to  solve  the  Xavier-Stokes 
equations  numerically. 

The  majority  of  the  published  numerical  studies  for  intermediate  Reynolds 
numbers  have  focused  on  flows  past  a  single  particle  and  are  thus  relevant  only  at  low 
particle  concentration.  In  regions  of  large  concentration,  the  drag  coefficient  is 
significantly  different  from  thac  of  an  isolated  particle  at  the  same  Reynolds  number, 
and  the  lift  and  moment  (torque)  coefficients  have  finite  values.  In  order  to 
understand  the  behaviour  of  a  particle  in  a  large-concentration  region,  studies  of  the 
interactions  amongst  particles  are  required.  Unfortunately,  in  practice,  the  spacial 
arrangement  of  particles  or  droplets  in  a  region  of  large  concentration  is  quite 
complex  and  subject  to  uncertainty,  and  calculations  involving  the  entire  region  are 
at  present  not  feasible.  In  order  to  develop  statistical  approaches,  information  about 
individual  droplet  or  particle  interactions  is  needed. 

The  study  of  droplet  or  particle  arrays,  particularly  in  flows  at  intermediate 
Reynolds  number,  is  relatively  new  (Patnaik  1986).  A  particle  array  as  discussed  by 
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Sirignano  (1983)  consists  of  a  few  particles  in  a  well-defined  geometrical  arrangement 
or  a  large  number  of  particles  in  a  periodic  configuration.  These  arrays,  although 
artificial,  can  provide  information  on  particle-particle  interactions  and  their  effects 
on  the  ambient  conditions  in  the  vicinity  of  the  particle.  The  simplification  in  the 
geometry  allows  a  detailed  and  rigorous  analysis. 

Tal,  Lee  &  Sirignan?!  (1983)  studied  the  hydrodynamics  and  heat  transfer  in 
assemblages  of  solid  spheres  in  a  steady  flow  at  Reynolds  number  100  (based  on  the 
particle  diameter  and  relative  velocity).  Their  method  took  advantage  of  the  periodic 
nature  of  an  infinite  array  of  spheres.  They  considered  several  spheres  in  tandem,  and 
foimd  a  trend  of  decreasing  drag  coefficients  in  the  streamwise  direction.  Tal,  Lee  & 
Sirignan^  (1984)  also  studied  the  interaction  of  two  sohd  spheres  in  tandem  in  a 
steady  uniform  fiow  at  =  40  for  two  different  spacings  using  bispherical 
coordinates  and  indicated  that  the  drag  coefficient  of  either  sphere  is  always  less  than 
that  of  a  single  isolated  sphere  and  that  the  reduction  is  much  greater  for  the 
downstream  sphere. 

Patnaik  (1986)  investigated  the  interaction  of  two  vaporizing  droplets  in  tandem 
at  Re  =  50  and  100  for  Lnterdroplet  spacing  equal  to  4.25  diameters  using  the 
doMmstream  solution  of  the  lead  droplet  as  the  inflow  conditions  for  the  downstream 
droplet.  He  found  that  for  both  Reynolds  numbers,  the  drag  coefficient  of  the  trailing 
droplet  is  lower  than  that  of  the  leading  droplet. 

Raju  &  Sirignano  ( 1990)  studied  the  interactions  between  two  moving  vaporizing 
droplets  in  tandem  at  50  ^  ^  200  for  a  range  of  spacings  and  droplet  radii  ratios. 

They  found  that  the  drag  coefficients  for  both  droplets  are  less  than  that  of  a  solid 
sphere  and,  for  the  same  Reynolds  number,  the  trailing  droplet  has  lower  drag. 
Chiang  &  Sirignano  (1991  a,  6)  ertended  this  study  to  include  variable  properties  and 
two  and  three  droplets.  With  three  droplets,  the  difference  in  drag  coefficients 
between  the  second  and  third  droplets  is  much  smaller  than  for  the  first  and  second 
droplets. 

All  of  the  above  studies  employed  axisymmetric  calcxiiations.  Recently,  some 
numerical  studies  have  been  performed  for  three-dimensional  flows  over  a  single  solid 
sphere.  Dandy  &  Dwyer  (1989)  obtained  three-dimensional  numerical  solutions  for 
steady,  uniform  shear  flow  past  a  fixed,  heated  spherical  particle  over  a  range  of 
Reynolds  numbers  (0.1  ^  Re  <  100)  and  dimensionless  shear  rates  (0.005  ^  a  ^  0.4). 
They  found  that  at  a  fixed  shear  rate,  the  lift  coefficient  is  approximately  constant 
over  a  wide  range  of  intermediate  Reynolds  numbers,  and  the  drag  coefficient  also 
remains  constant  when  normalized  by  the  drag  for  a  sphere  in  imiform  flow. 
Tomboulides,  Orszag  &  Kamiadakis  (1991)  performed  a  numerical  study  of  three- 
dimensional  flow  past  a  sphere  using  a  spectral  element  method  for  30  ^  Re  ^  1000 
and  discussed  steady  axisymmetric  states  and  unsteady  states  with  three- 
dimensional  vortex  shedding. 

Three-dimensional  flow  interactions  between  droplets  or  particles  at  finite 
Reynolds  number  have  not  yet  been  studied.  As  a  first  step  towards  understanding 
the  three-dimensional  interactions  in  large  concentration  of  particles,  we  investigate 
flow  interactions  between  two  identical  (solid  or  liquid)  spheres  which  are  held  fixed 
side  by  side  against  the  uniform  stream  at  Reynolds  numbers  50,  100,  and  150.  We 
determine  the  effects  of  three-dimensional  interactions  on  the  lift,  moment,  and  drag 
coefficients  as  a  function  of  the  dimensionless  distance  between  the  two  spheres  and 
Reynolds  number.  Some  novel  phenomena  in  the  near  wake  are  discovered  as  the  gap 
between  the  two  spheres  decreases.  Our  results  are  also  compared  with  available 
experimental  and  numerical  data. 


Three-dimensional  flow  over  two  spheres 
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Figttrz  1.  Flow  geometry  and  coordinates. 


2.  Problem  statement  and  numerical  solution 

We  consider  a  steady  three-dimensional  incompressible  laminar  flow  of  a 
Newtonian  fluid  past  two  identical  (solid  or  liqtiid)  spheres  held  fixed,  with  the  line 
connecting  the  sphere  centres  normal  to  a  uniform  stream,  as  shown  in  figure  1 ; 
denotes  the  distance,  normalized  by  the  sphere  radius,  from  the  sphere  centre  to  the 
r,  y  symmetry  plane  between  the  two  spheres.  Far  upstream,  the  flow  is  uniform  with 
constant  velocity  i  parallel  to  the  r-axis.  Two  symmetry  planes  are  noted  in  figure 
1 ;  the  (x,  3)-plane  containing  the  centres  of  the  two  spheres  and  the  (z,  y)-plane 
located  at  z  =  —  dj  midway  between  the  sphere  centres. 

We  note  that  asymmetry  of  the  flow  past  neighbouring  bluff  bodies  might  occur 
at  lower  Reynolds  number  than  that  of  the  first  temporal  instability  for  a  single 
body.  An  example  is  the  flow  past  two  cylinders  where  flow  asymmetry  would  occur 
well  before  the  first  temporal  instability  for  a  single  cylinder.  In  contrast,  we  expect 
that,  for  the  flow  past  two  spheres,  flow  asymmetry  would  occur  nearly’ 
simultaneously  with  the  temporal  instability  since  the  interaction  between  the  two 
halves  (each  containing  one  sphere)  of  the  flow  field  should  be  stronger  for  cylinders 
than  for  spheres.  That  is,  the  flow  between  the  two  spheres  is  less  constrained  than 
the  flow  between  the  two  side-by-side  cylinders. 

Two  coordinate  systems  are  used  in  our  formulation :  the  Cartesian  coordinates  (x, 
y,  z)  and  the  non-orthogonal  generalized  coordinates  (^,  y,  Q.  The  origin  of  the  former 
coincides  with  the  sphere  centre;  and  $  is  the  radial,  y  the  angxilar,  and  ^  the 
azimuthal  coordinate.  Owing  to  symmetry,  the  physical  domain  is  reduced  to  one 
quarter  of  an  ellipsoid-like  space. 

The  domain  of  the  external  flow  is  bounded  by  1  $  ^  ^  .Vp  1  ^  7  <  X,  1  ^  ^  ^ 
where  i  =  I  and  correspond,  respectively,  to  the  sphere  surface  and  the  far-field 
boundary  surrounding  the  sphere ;  7  =  1  and  X  denote,  respectively,  the  positive  z- 
axis  and  the  negative  z-axis;  ^  =  I  and*Vj  refer,  respectively,  to  the  (x,  z)-plane  in  the 
positive  r-direction  and  the  {i.z)-plane  in  the  negative  r-direction. 

The  domain  of  the  internal  flow  is  1  ^  ^  1  ^  7,  <  Nj,  1  ^  ^  *Vj.  =  1  and 

iVij  correspond  to  the  centre  and  the  surface  of  the  sphere,  respectively.  7;  =  1  and 
X  denote  the  positive  z-axis  and  the  negative  z-axis,  respectively.  =  1  and  N,  refer 
to  the  (z,  z)-plane  in  the  positive  z-direction  and  the  (z,z)-plane  in  the  negative 
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z-direction,  respectively.  Within  the  liquid  sphere,  =  constant  are  a  family 
of  concentric  spherical  surfaces.  Uniform  spacing  =  Stj  =  =  1)  is  used,  for 

convenience,  for  the  generalized  coordinates  in  both  flows. 

The  non-orthogonal  generalized  coordinate  s^rstem  of  the  external  flow  can  be 
easily  adapted  to  three-dimensional  arbitrary  geometries.  We  solve  the  continuity 
equation  and  the  time-dependent  Navier-Stokes  equations  and  relax  them  to  the 
steady-state  solution,  as  will  be  discussed  in  detail  in  §2.2. 


2.1.  Governing  eqvxiiicms  and  boundary  conditions 

Since  one  of  our  goals  is  to  study  the  flow  interaction  with  liquid  spheres,  we  present 
the  equations  for  the  flows  inside  and  outside  the  spheres.  However,  for  flow 
interactions  with  solid  spheres,  only  the  external  flow  equations  are  solved.  The 
continuity  and  momentum  equations  inside  and  outside  a  sphere  and  the  boundary 
conditions  are  non-dimensionalized  using  the  sphere  radius  a#  as  the  characteristic 
length  and  as  the  characteristic  velocity; 

external  flow 


V-U  =  0, 


DU 

Dt 


(1) 

(2) 


internal  flow 


V.P^  =  0, 


Dt 


(3) 

(4) 


The  governing  equations  are  written  with  respect  to  the  generalized  coordinates 
(^1 V'  0>  which  allows  a  three-dimensional  body  of  arbitrary  shape  to  be  treated.  The 
numerical  integration  is  performed  using  a  cubic  computational  mesh  with  equal 
spacing  (S^  =  Sy  =  1)  (Anderson,  Tannehill  &  Fletcher  1984). 

The  conditions  at  the  interface,  ^  =  1  or  are  derived  by  requiring 

continuity  of  the  shear  stresses  and  tangential  velocities.  Because  no  fluid  is  allowed 
to  cross  the  surface  of  the  liquid  sphere,  the  normal  velocities  at  the  interface  are  zero 
in  both  flows.  An  interface  condition  for  the  pressure  is  obtained  from  the  momentum 
equation.  Since  the  interface  is  always  spherical,  it  is  more  convenient  to  cast  these 
conditions  in  terms  of  spherical  coordinates  (r,  $,  (p) : 


+— 


“^re.s  — 

(5) 

(6) 

=  Vi.s.s. 

(7) 

(8) 

K..=  Vi.r.,=0, 

(9) 

cr  \  Sr  y  S0  sm  0  C0  _ 

(10) 

where  the  subscript  5  denotes  the  surface  of  the  liquid  sphere.  For  the  solid-sphere 
case,  the  no-slip  condition  is  enforced  on  the  sphere  surface. 


Three-dimensicmal  jiow  over  two  spheres 
The  external  flow  boundary  conditions  are ; 


p  =  0,  u=l,  y  =  0,  iu  =  0  at  f  =  .V^  except  at  z  = —do, 


^p  _^u  _  dv 
5^“ 


0.  u;  =  0  at  z 


-d 


o> 


(11) 

(12) 


dv  fiit  cw 

•^  =  T-  =  — =  0,  y  =  0  at  7  =  1  and  N,,  (13) 

07  27  C7 

5v  6u  cw 

^  =  ^  =  •^  =  0,  y  =  0  at  ^=1  and  N,,  (14) 

where  u,  v,  and  w  are  the  velocities  of  the  external  flow  in  the  x-,  y-,  and  s-directions, 
respectively,  p  is  the  pressure,  and  the  subscript  I  denotes  the  internal  flow. 
Equations  (12)  and  (14)  correspond,  respectively,  to  the  symmetry  conditions  in  the 
X,  y  symmetry  plane  between  the  spheres  and  the  x,  z  symmetry  plane  containing 
the  centres  of  the  spheres.  Equation  (13)  expresses  the  no-flux  boimdary  condition 
for  p,  u,  and  w  on  the  axes  7  =  1  and  7  =  -V,.  and  zero  v-velocity  in  the  x,  z  symmetry 
plane  containing  the  two  axes. 

The  internal  flow  boundary  conditions  are 


where  (15)  and  (16)  correspond,  respectively,  to  the  no-flux  boundary  conditions  at 
the  centre  of  the  droplet  and  on  the  axes  7,  =  1  and  iV„  and  zero  u, -velocity  in  the* 
X,  2  symmetry  plane  containing  the  two  axes.  Equation  (17)  prescribes  the  symmetry 
condition  on  the  x,  2  symmetry  plane  containing  the  centres  of  the  spheres. 

The  dimensional  drag,  lift,  and  moment  coefiBcients  are  evaluated  from 


Eo  =  I  — pn-i(iS-l- 1  n-r-idJS, 

(18) 

Js  Js 

—  {  —pn-kdS+  n-t-kd3, 

(19) 

Js  Js 

Jf  =  J  rxrdS, 

(20) 

where  S  denotes  the  surface  of  the  sphere,  n  is  the  outward  unit  normal  vector  at  the 
surface,  r  is  the  position  vector  from  the  centre  of  the  sphere,  and  t  is  the  viscous 
stress  tensor.  Equations  (18)-(20)  are  evaluated  on  the  side  of  the  external  flow.  The 
lift  force  is  assumed  positive  when  it  is  directed  toward  the  positive  2-axis.  Owing  to 
symmetry,  only  the  y-componen*:  of  the  moment  is  non-zero  and  is  assumed  positive 
in  the  counter-clockwise  direction. 
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Non-dimensional  coefficients  of  drag,  lift,  and  moment  are  defined  respectively  as 


Co  = 

(21) 

Ci=FJ{y)U\,Kal), 

(22) 

(23) 

2.2.  Numerical  solution 

We  have  developed  a  three-dimensional  implicit  finite-difference  algorithm  to  solve 
simultaneously  the  set  of  discretized  partial  differential  equations.  The  method  is 
based  on  an  Altemating-Direction-Predictor-Corrector  (AD PC)  scheme  to  solve  the 
time-dependent  Navier-Stokes  equations.  ADPC  is  a  slight  variation  of  Altemating- 
Direction-Implicit  (ADI)  method.  It  is  first-order  accurate  in  time  but  is  effective 
and  implemented  easily  when  embedded  in  a  large  iteration  scheme  (Patnaik  1986). 
The  control  volume  formulation  is  used  to  develop  the  finite-difference  equations 
from  the  governing  equations  with  respect  to  the  generalized  coordinates 
One  of  the  advantages  of  the  control  volume  formulation  is  that  all  dependent 
variables  are  conserved  over  a  single  control  volume,  and  hence  the  whole  domain 
regardless  of  the  grid  fineness.  An  important  part  of  solving  the  Navier-Stokes 
equations  in  primitive  variables  is  the  calculation  of  the  pressure  field.  In  the  present 
work,  a  pressure  correction  equation  is  employed  to  satisfy  indirectly  the  continuity 
equation  (Anderson  et  al.  1984).  The  pressure  correction  equation  is  of  the  Poisson 
type  and  is  solved  by  the  Successive-Over-Relaxation  (SOR)  method. 

The  overall  solution  procedure  is  based  on  a  cyclic  series  of  guess-and-correct 
operations.  The  velocity  components  are  first  calculated  from  the  momentum 
equations  using  the  ADPC  method,  where  the  pressure  field  at  the  previous  time  step 
is  employed.  This  estimate  improves  as  the  overall  iteration  continues.  The  pressiire 
correction  is  calculated  from  the  pressure  correction  equation  using  the  SOR  method, 
and  new  estimates  for  pressure  and  velocities  are  obtained.  This  process  continues 
until  the  solution  converges  at  each  time  step.  Por  the  flow  past  liquid  spheres,  the 
same  procedure  is  performed  in  the  flow  inside  the  sphere.  The  governing  equations 
of  motion  in  each  flow  are  solved  in  an  interactive  sequence  through  the  interface 
boxmdary  conditions  until  convergence  is  achieved  for  each  time  step  of  the 
calculation. 

The  generation  of  the  computational  grid  for  the  external  multi-sphere  flows  is  an 
essential  part  of  the  solution  procedure.  We  generate  the  three-dimensional  grid 
efficiently  by  choosing  the  outer  boundary  of  the  physical  domain  to  be  axisymmetric 
about  the  line  connecting  the  centres  of  the  two  spheres  and  constructing  an 
axisymmetric  grid.  The  axisymmetric  grid  is  generated  on  the  (^,i7)-plane  including 
the  line  that  connects  the  centres  of  the  two  spheres  and  by  using  a  hybrid  method 
of  algebraic  and  differential  equation  methods.  First,  we  choose  the  outer  boundary 
of  the  computational  grid  as  the  x.  y  symmetry  plane  and  an  incomplete  ellipse  whose 
centre  is  located  at  the  centre  of  the  sphere,  as  shown  in  figure  1.  We  then  generate 
a  family  of  quarter-ellipses  on  the  right  side  of  the  domain,  another  family  of  quarter- 
ellipses  on  the  left  side,  and  also  a  family  of  straight  lines  emanating  from  the  centre 
of  the  sphere.  Grid  density  is  controlled  by  the  stretching  function  developed  by 
Vinokur  (1983).  This  is  followed  by  solving  the  quasi-linear  elliptic  system  of 
differential  equations  using  the  SOR  method  with  a  few  (not  more  than  two) 
iterations  in  order  to  smooth  the  grid  lines  generated  by  the  algebraic  method.  In 
figures  2  (a-c),  we  present  a  cross-section  of  the  three-dimensional  grid  at  d,,  =  21,  10, 
and  3. 
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FiauBK  2.  Cross-section  of  the  three-dimensional  grid  system  for  (a)  i,  =  21;  [b)  dj  =  10; 

(c)  ^0  ~  3- 


xVjX.VjX.Vj 

Cgy 

Cl 

Re  =  50 

20x21  X  21 

0.706 

0.951 

1.657 

30x31  X  31 

0.683 

0.934 

1.617 

40x41x41 

0.676 

0.929 

1.605 

1.58 

It 
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20x21x21 

0.553 

0.593 

.  1.148 

30x31  X  31 

0.532 

0.582 

1.114 

40x41x41 

0.524 

0.581 

1.105 

1.09 

Table  1.  Drag  coefficients  as  a  function  of  grid  density  at  Re  =  30  and  100  where  *  denotes  the 
data  from  Roos  &  Willmarth  (1971)  and  also  Clift  et  al.  (1978). 


3.  Results  and  discussion 

In  §3.1,  we  test  the  accuracy  of  the  full  three-dimensional  solution  procedure  by 
predicting  the  axisymmetric  flow  over  a  single  (solid  and  liquid)  sphere.  In  §3.2,  we 
discuss  the  three-dimensional  interactions  between  two  solid  spheres,  and  in  §3.3,  we 
examine  the  three-dimensional  interactions  between  two  liquid  spheres. 

3.1.  Flow  over  a  -si-ngle  sphere 

We  discuss  the  flow  generated  by  an  impulsively  started  solid  sphere  in  a  quiescent 
fluid  at  two  Reynolds  numbers :  50  and  100.  The  time-dependent  solution  converges 
as3rmptotically  to  a  steady  state,  which  is  in  excellent  agreement  with  available 
experimental  data  and  correlations  ais  shown  in  tables  1  and  2.  Table  1  lists  the  drag 
coeflScients  as  a  function  of  the  computational  grid  density  at  Reynolds  numbers  50 
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Figitbe  4.  Velocity  vector  fields  of  (a)  external  flow;  (6)  internal  flow  for  the  axisymmetric  flow 

past  a  liquid  sphere  at  Re  =  100. 


The  axisymmetric  test-run  for  a  solid  sphere  at  Reynolds  number  100  with  the 
40  X  41  X  41  grid  required  a  dimensionless  time  step  of  Ai  =  0.002  and  a  total  time  of 
3.75  Cray  Y-MP/8-64  cpu  hours.  Each  time  step  takes  about  1.8  cpu  seconds. 

We  also  solved  the  flow  generated  by  an  impulsively  started  liquid  sphere  in  a 
quiescent  fluid,  and  the  results  were  in  good  agreement  with  available  numerical 
studies.  Figure  3  (a,  b)  shows  the  steady -state  velocity  vector  distributions  of  the  flow 
past  a  liquid  sphere  with  viscosity  ratio  of  25  and  density  ratio  of  300  (internal  to 
external  fluid)  at  Reynolds  number  30.  We  observe  a  closed-streamline  wake 
detached  from  the  liquid  sphere,  and  thus  no  secondary  recirculating  flow  is  found 
in  the  liquid  phase  (Rivkind  &  Ryskin  1976  ;  Clift  et  al.  1978).  In  figure  4  (a,  b)  we 
show  the  velocity  vector  distributions  for  the  same  parameters  as  above  except  that 
the  Reynolds  number  is  100.  It  is  interesting  to  note  that  a  second  circulatory  flow 
develops  in  the  liquid-sphere  stem  region.  This  behaviour  was  observed  in  an  earlier 
study  by  Rivkind  &  Ryskin  (1976)  where  a  stream  function-vorticity  formulation 
was  employed.  Rivkind  &  Ryskin  (1976)  indicated  that  when  the  density  ratio  is 
much  greater  than  the  viscosity  ratio,  i.e.  the  Reynolds  number  inside  the  liquid 
sphere  is  much  greater  than  the  Reynolds  number  outside,  a  second  circulatory  flow 
possibly  occurs  in  the  liquid-sphere  stem  region.  The  axisymmetric  test  run  for  the 
liquid  sphere  required  half  the  time  step  for  the  solid  sphere  and  was  about  3.4  times 
slower  because  of  the  numerical  interaction  between  the  internal  and  external  flows. 

For  the  interactions  between  two  spheres,  an  ellipsoid-Uke  domain  is  chosen  in 
order  to  take  into  account  the  interactions  when  the  two  spheres  are  far  away  from  ’ 
each  other.  As  shown  in  figure  1,  a  longer  outer  boundary  =  25  is  chosen  in  the  :- 
direction,  and  =  21  is  chosen  in  the  o'-direction,  where  cr  =  (i*-t-y-)V  The  results 
using  the  above  ellipsoid  outer  boundary  for  a  single  sphere  were  identical  in  the 
steady  axisymmetric  flow  calculations  to  those  using  the  spherical  outer  boundary. 

The  results  of  the  30  x  31  x  31  grid  differ  by  only  0.8%  in  the  drag  coefficients  and 
0.4%  in  the  separation  angles  from  those  of  the  40  x  41  x  41  grid  as  shown  in  tables 
1  and  2.  The  percentage  difference  was  calculated  ais  follows.  Let  the  result  of  the 
30x31x31  grid  be  and  the  result  of  the  40x41x41  grid  be  5,.  Then,  the 
percentage  difference  is  {S^  —  S^)/S^.  Thus,  we  chose  the  medium-size  grid 
30x31x31,  and  15x31x31  inside  the  liquid  sphere,  for  the  following  of 
calculations.  A  typical  run  for  the  solid  sphere  with  the  30  x  31  x  31  grid  required  0.8 
cpu  hours  on  the  Cray  Y-MP/8-64;  the  liquid  sphere  run  with  the  same  grid  required 
2.7  cpu  hours.  15  runs  were  made  for  each  Remolds  number  to  cover  the  range  of 
1.5  $  dg  <  25. 

We  did  not  perform  calculations  for  Reynolds  number  higher  than  150.  because  it 
is  known  that  the  wake  becomes  unstable  at  a  Reynolds  number  in  the  range  of 
130-220  for  a  single  sphere  (Taneda  1956;  Goldburg  &  Florsheim  1966  ;  Roos  & 
Willmarth  1971 ;  Nakamura  1976;  Kim  A  Pearlstein  1990)  and  our  present  goal  is  to 
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FiGtTBB  5.  Sketch  of  typical  streamlines  over  one  of  the  two  solid  spheres  in  the  principal  plane 
(z, /-plane)  at  Re  =  100  for  (a)  axisymmetric  flow;  (6)  d,  =  2;  (c)  d,  =  1.5. 

obtain  steady-state  solutions.  For  the  solution  of  a  flow  including  the  three- 
dimensional  unsteady  wake,  a  complete  computational  domain  (i.e.  encompassing 
the  two  spheres  without  a  sjTnmetry  plane)  and  periodic  boundary  conditions  in 
^-direction  will  be  necessary. 

3.2.  Interactions  of  two  solid  spheres 

3.2.1.  Flow  structure 

In  order  to  illustrate  better  the  fluid  motion,  we  consider  the  flow  field  in  the  (z:,  z)- 
plane  of  symmetry,  which  is  defined  as  the  principal  plane,  where  the  narrowest 
path  between  the  two  spheres  occurs,  hence  the  strongest  interactions  between  them. 

Figure  5(a)  displays  a  sketch  of  the  actual  streamlines  around  a  single  sphere  in  the 
principal  plane  at  Reynolds  number  100.  As  expected,  two  identical  counter-rotating 
vortices  (corresponding  to  a  single  vortex  ring)  exist  in  the  wake,  and  the 
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downstream  stagnation  point  is  located  on  the  axis  of  symmetry.  Line  I-II 
connecting  the  front  and  rear  stagnation  points  in  the  standard  axisymmetric  dow 
over  a  single  sphere  will  be  used  as  a  reference  line :  we  refer  to  the  region  above  this 
line  as  the  ‘  top  ’  or  '  upper  ’  region  and  that  below  as  the  ‘  bottom '  or  ‘  lower  ’  region. 

Figure  5(6)  displays  a  sketch  of  the  actual  streamlines  around  one  of  the  two 
spheres  in  the  principal  plane  at  Reynolds  number  100.  The  two  spheres  are 
separated  by  a  distance  dj  =  2.  Owing  to  the  blockage  of  the  flow  in  the  gap  between 
the  two  spheres,  the  streamlines  diverge  away  from  the  x,  y  symmetry  plane  (located 
at  z  =  —d^)  as  they  approach  the  front  stagnation  region.  Thus,  the  stagnation 
streamline  of  the  single-sphere  case  (I-Sj  in  figure  5a)  no  longer  intersects  the  sphere, 
and  another  streamline  closer  to  the  symmetry  plane  meets  the  sphere  to  form  the 
new  front  stagnation  streamline  at  point  P.  As  a  consequence,  the  fluid  particles 
move  faster  in  the  lower  left  region  around  the  sphere  than  in  the  upper  left  region, 
and  tins  causes  the  pressure  in  the  lower  left  region  to  be  lower  than  that  in  the  upper 
left  region.  The  resulting  pressure  difference  between  the  upper  and  lower  left  regions 
is  higher  than  that  between  the  bottom  of  the  sphere  and  the  narrow  path.  This 
pressure  imbalance,  which  will  be  discussed  in  §3.2.2,  causes  repulsion  of  the  two 
spheres.  The  contribution  of  shear  stress  differences  to  the  repulsion  will  also  be 
discussed  in  §3.2.2. 

Figure  5(6)  shows  an  interesting  streamline  pattern  in  the  wake  region.  Two 
counter-rotating  eddies  exist  in  the  wake  but  their  configuration  is  quite  different 
from  that  for  axisymmetric  flow.  The  lower  eddy  is  formed  by  the  fluid  separating 
on  the  lower  portion  of  the  sphere  as  in  the  case  of  axisymmetric  flow.  The  upper 
eddy  is  not  formed  by  the  fluid  separating  on  the  upper  portion  of  the  sphere,  but 
rather  by  the  fluid  turning  around  the  lower  eddy  and  being  entrained  by  the  upper 
flow.  This  upper  eddy  is  detached  from  the  sphere.  A  portion  of  the  fluid  moving 
around  the  bottom  of  the  sphere  passes  between  the  detached  upper  eddy  and  the 
sphere.  The  streamline  A-B  encompassing  the  upper  eddy  intersects  itself,  and  the 
intersection  point,  C,  designated  as  the  downstream  stagnation  point,  is  shifted, 
toward  the  x.y  symmetry  plane.  Both  eddies  are  smaUer  than  those  of  the 
axisymmetric  flow.  These  new  features  can  be  explained  as  follows.  The  pressure 
above  the  upper  wake  is  less  than  that  below  the  lower  wake  owing  to  the  increased 
acceleration  of  the  fluid  in  the  narrow  path  between  the  two  spheres  (as  will  be  shown 
in  figure  7).  Thus,  the  flxiid  particles  turning  aroimd  the  lower  eddy  are  pushed  into 
the  upper  region  of  the  wake.  The  pressure  distribution  around  the  sphere  will  be 
discussed  further  in  §3.2.2. 

Figure  5  (c)  shows  a  sketch  of  the  actual  streamline  pattern  at  Reynolds  number 
100  for  the  case  d,  =  1.5.  The  shifting  of  the  front  stagnation  streamline  and 
stagnation  point  toward  the  x,  y  symmetry  plane  are  more  obvious  here  than  in  the 
previous  case  of  d,  =  2.  The  significant  difference  is  in  the  wake  region  where  both  the 
upper  eddy  and  downstream  stagnation  point  vanish.  Fluid  particles  separating  on 
the  upper  portion  of  the  sphere  move  downstream  without  returning  (streamline 
r>-E).  On  the  other  hand,  fluid  particles  turning  around  the  lower  eddy  move  into  the 
upper  region  of  the  wake  until  they  reach'near  the  upper  separation  point,  D.  and 
then  move  downstream  in  an  S-shaped  path  (streamline  F-G)  without  returning  to 
form  an  eddy.  The  lower  eddy  shrinks  as  the  two  spheres  become  closer,  and  the 
pressure  difference  between  the  top  and  bottom  of  the  wake  is  larger. 

It  is  interesting  to  examine  the  changes  in  the  separation  region  at  the  sphere 
surface  for  the  cases  d,  =  2  and  1.5.  More  specifically,  we  examine  the  behaviour  of 
the  circle  of  intersection  of  the  wake  and  the  sphere.  Our  results  show  that  the  circle 
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Figctre  6.  Velocity  vector  fields  over  one  of  the  two  solid  spheres  in  the  principal  plane 
corresponding  to  (a)  figure  5(6),  =  2;  (6)  figure  5(c),  d,  =  1.5. 


is  slightly  shifted  toward  the  x,  y  symmetry  plane  due  to  the  decreased  pressure  in  the 
gap  region  with  respect  to  that  in  the  wake  lower  region.  This  shifting  produces 
separation  angles  at  the  top,  middle,  and  bottom  of  the  sphere  with  values  of  123.1®, 
126.5®,  and  126.2®  respectively  for  the  case  =  1.5,  where  the  angles  are  measured 
from  the  front  stagnation  point  of  the  axisynametric  flow  case.  This  is  in  contrast 
with  an  angle  of  125.7®  at  all  positions  for  a  single  sphere. 

Figure  6(a)  shows  the  computed  velocity  vector  field  corresponding  to  figure  5(6), 
dj  =  2.  The  velocity  vectors  upstream  of  the  front  stagnation  point,  S^,  for 
axis3mimetric  flow  point  downwards  away  from  the  x,y  symmetry  plane.  This 
indicates  that  the  front  stagnation  point  (P  in  figure  56),  is  shifted  toward  the  x,  y 
symmetry  plane.  Similar  behaviour  occurs  at  the  rear  stagnation  point  (Q  in  figure 
56).  Figure  6(6)  shows  the  computed  velocity  vector  field  corresponding  to  figure 
5(c),  do  =  1.5. 

One  of  the  advantages  of  the  velocity  vector  plot  is  that  it  shows  clearly  the' 
relative  magnitude  of  velocity  in  the  flow  field,  e.g.  the  smaller  velocity  in  the  wake 
region  compared  to  that  outside  the  wake  is  seen  in  figure  6  (a,  6).  However,  it  is 
difficult  to  obtain  streamlines  directly  from  the  tangents  of  the  velocity  vector  plot. 

A  stream  function  ^  cannot  be  defined  and  calculated  from  the  velocity  in  the 
principal  plane  owing  to  the  existence  of  a  divergence  associated  with  the  third 
component  of  velocity.  Nevertheless,  for  descriptive  purposes  only,  it  is  convenient 
to  use  the  algorithm 

Vi^ps(r. <9o)  =  +  f  -"“odr  (24) 

■'’■o 

to  present  appro.ximations  to  the  streamline  pattern.  It  is  understood  that  since  a 
true  stream  function  does  not  exist,  the  pseudo-stream  function  is  dependent  upon 
the  integration  path.  The  above  algorithm  specifically  involves  only  radial 
integration ;  can  be  recovered  by  differentiation  of  this  function,  but  u,.  cannot  be 
recovered.  The  streamlines  presented  in  figures  o(a-c)  were  based  on  this  algorithm. 

Phenomena  in  the  wake  similar  to  those  described  above  have  been  found  in  a  few 
previous  studies.  Rosfjord  (1974)  obtained  results  similar  to  those  in  figures  5  (6)  and 
5(c)  in  his  experimental  and  numerical  studies  of  the  recirculating  flow  region 
between  two-dimensional  parallel  separated  jets.  He  found  that  for  velocity  ratios 
between  two  jets  equal  to  1.11  and  1.25  (upper  to  lower),  two  eddies  exist  near  the 
injector  face,  but  the  upper  eddy  is  detached  from  the  injector  face,  and  a  portion  of 
the  fluid  originating  at  the  lower  jet  is  entrained  by  the  upper  jet,  passing  between 
the  detached  upper  eddy  and  the  injector  face.  He  also  reported  that  for  a  velocity 
ratio  of  1.4,  only  the  lower  eddy  existed  and  a  complete  entrainment  of  the  weaker 
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Fioitb*  7.  (a)  Pathlines  of  two  fluid  particles  (A  and  B)  slightly  above  the  principal  plane 
(y  =>  0.001)  of  the  two  sobd  spheres.  (6,c)  Pathline  of  a  fluid  particle  (C)  whose  initial  position  was 
(*,  y,z=*  0.996, 0.575,  0) :  (6)  a  view  from  top  of  the  sphere  looking  toward  the  z,  y  symmetry  plane , 
(c)  a  side  view  looking  normal  to  the  z,  y  symmetry  plane. 
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Figtoe  8.  Distribution  of  the  pressure  coefBcient  around  the  solid  sphere  in  the  principal  plane 

at  Ri  =  100  for  i,  =  1.5. 

jet  was  observed.  In  particular,  the  flow  S-shaped  loop  near  the  stronger  jet  was 
clearly  indicated.  Recently,  Dandy  k  Dwyer  (1990)  also  found  a  flow  pattern  similar 
to  figure  6  (6)  in  their  numerical  study  of  steady  uniform  shear  flow  past  a  single  solid 
sphere. 

In  order  to  facilitate  the  visualization  of  the  three-dimensional  character  of  the 
flow  in  the  wake  region  discussed  in  detail  in  §3.2.1,  we  present  the  pathlines  of 
selected  fluid  particles  in  figure  7  (a-c),  at  Reynolds  number  100  for  =  2,  where  the. 
free-stream  direction  is  from  left  to  right.  The  pathiines  x(z^,y^,ZQ,t),  where  the 
subscript  0  denotes  initial  particle  location,  were  obtained  by  solving  three  coupled 
ordinary  differential  equations  dx/di  =  u{x),  via  a  fourth-order  Runge-Kutta 
method. 

We  first  selected  two  fluid  particles  (A  and  B)  slightly  above  the  principal  plane 
(y  =  0.001)  separated  by  a  small  distance  (much  smaller  than  the  sphere  radius) 
in  the  wake  region.  Figure  7  (a)  shows  that  particle  A  follows  an  S-shaped  pathline. 
On  the  other  hand,  particle  B  follows  a  closed-loop  pathline  as  was  discussed  in  the 
previous  section. 

We  then  examined  the  pathline  of  a  fluid  particle  C  whose  initial  position  (ar,,, 

Zg  =  0.995,  0.575,  0)  was  in  the  wake  region  but  above  the  principal  plane.  Figure  7  (6) 
(a  view  from  top  of  the  sphere  looking  toward  the  x,  y  symmetry  plane)  shows  that 
the  fluid  particle  C  first  follows  a  helical  pathline  as  it  approaches  the  x,  y  symmetry 
plane  and  then  moves  downstream.  Figme  7  (c)  is  a  side  view  (looking  normal  to  the 
X,  y  symmetry  plane)  of  that  pathline. 

3.2.2.  Pressure  and  shear  stress  distribution 

Figure  8  shows  the  pressure  coefficient,  -{p—Pja)/pC%>  around  one  of  the  spheres 
in  the  principal  plane  at  Reynolds  number  100  ford,  =  1.5.  On  average,  the  pressure 
is  higher  on  the  top,  contributing  to  a  positive  lift  force.  The  pressure  on  the  bottom 
of  the  sphere  is  lower  between  0®  and  87.6°  and  also  slightly  lower  between  159°  and 
180°  in  the  wake  region  than  that  on  the  top  of  the  sphere,  where  the  angle  is 
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Figttbe  9.  Distribution  of  shear  stress  coefficient  around  the  solid  sphere  in  the  principal  plane 

at  Re  =  too  for  dj  =  1.5. 
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Ftotmi  10.  Tangential  velocity  profile  on  the  bottom  and  top  of  the  solid  sphere  in  the 
principal  plane  at  (a)  d  =  36°;  (6)  d  =  84°; - ,  single  sphere; . ,  bottom; - ,  top. 

measured  from  the  front  stagnation  point  [d  =  0)  of  the  axisymmetric  flow  case.  On 
the  bottom  of  the  sphere,  the  minimum  pressure  occurs  at  an  angle  less  than  90°.  On 
the  top  of  the  sphere,  the  minimum  pressure  occurs  at  an  angle  greater  than  90°  and 
is  lower  than  the  minimum  pressure  on  the  bottom  of  the  sphere.  The  maximum 
pressure  occurs  a  few  degrees  toward  the  symmetry  plane  measured  from  5  =  0. 
The  highest  pressure  in  the  wake  region  occurs  a  few  degrees  toward  the  z,y 
symmetry  plane  measured  from  d  =  k.  These  observations  indicate  that  the  front 
and  rear  stagnation  points  are  shifted  a  few  degrees  toward  the  x,  y  symmetry  plane. 

Figure  9  shows  the  tangential  shear  stress  coefBcient,  in  the  same  plane 

used  for  the  pressure  coefficient  in  figure  8.  Note  that  the  clockwise  direction  is 
considered  positive  for  the  shear  stress  on  the  top  of  the  sphere,  and  the 
counterclockwise  direction  is  considered  positive  for  the  shear  stress  on  the  bottom. 
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It  is  seen  that  the  shear  stress,  on  average,  is  higher  around  the  lower  part  of  the 
sphere  than  around  the  top.  In  particular,  the  magnitude  of  the  shear  stress  is  higher 
in  the  lower  regions  d  -  0  to  63.6  and  (9  =  165.5  to  it  than  on  the  top  of  the  sphere. 
It  is  also  important  to  note  that,  owing  to  their  inclinations  with  the  r-axis,  the  shear 
forces  on  these  two  lower  regions  contribute  to  both  the  lift  (parallel  to  the  ;-axis) 
and  drag  (parallel  to  the  z-axis),  whereas  the  shear  force  at  the  top  of  the  sphere 
contributes  mainly  to  the  drag.  Thus,  the  shear  forces  in  this  case  contribute,  along 
with  the  pressure  forces,  to  the  repulsion  of  the  two  spheres.  The  shear  stress  at 
6  =  0  is  not  zero  but  acts  counterclockwise,  and  the  shear  stress  at  ^  =  rc  is  not  zero 
but  acts  clockwise.  Therefore,  the  front  and  rear  stagnation  points  are  shifted  a 
few  degrees  toward  the  x,  y  symmetry  plane.  Another  interesting  feature  is  that  the 
separation  angle  where  the  shear  stress  vanishes  on  the  top  of  the  sphere  is  123.1°, 
but  that  angle  on  the  bottom  is  126.2°.  The  computed  separation  angle  of  125.7°  for 
the  axisymmetric  flow  case  with  the  medium-size  grid  was  shown  in  table  2.  Thus, 
the  reverse  flow  in  the  wake  region  is  shifted  upwards  toward  the  i,  y  symmetry 
plane. 

We  examine  next  the  tangential  velocity  profiles,  Ug[r),  at  two  different  fl-locations 
on  the  bottom  and  top  of  the  sphere  in  the  principal  plane  at  Reynolds  number  100 
for  dg  =  1.5.  Figure  10(a)  shows  the  tangential  velocity  profiles  aX  d  —  36°  on  the 
bottom  and  top  of  one  of  the  spheres,  in  addition  to  that  for  the  axisymmetric  case 
as  a  reference.  It  is  seen  that  the  maximum  velocity  at  the  bottom  is  higher  than  at 
the  top,  M  mentioned  earlier  in  the  discussion  of  figure  5(6).  It  is  also  seen  that  the 
velocity  gradient  at  the  sphere  surface  is  higher  at  the  bottom  than  at  the  top,  hence 
the  higher  shear  stress  at  the  bottom  as  explained  earlier  (see  figure  9).  Figure  10(6) 
shows  the  tangential  velocity  profiles  at  =  84°  on  the  bottom  and  top  of  the  sphere, 
in  addition  to  that  for  the  axisymmetric  case.  It  is  seen  that  now  the  velocity 
gradient  at  the  top  wall  is  26  %  higher  than  at  the  bottom  although  the  maximum 
velocity  at  the  top  is  only  3%  higher  than  the  maximum  value  at  the  bottom.  -The 
reason  is  that  the  boundary -layer  growth  at  the  top  is  limited  by  the  interaction  with 
the  boundary  layer  of  the  other  sphere. 

3.2.3.  Lift  coefficients 

In  the  following  discussion,  we  classify  the  proximity  of  the  two  spheres  into  three 
regimes;  close,  intermediate,  and /ar  separated,  depending  on  the  values  of  d^  and 
Reynolds  number. 

Figure  11  (a-c)  show  the  total  lift  coefficient  and  the  lift  coefficients  due  to  viscous 
and  pressure  contributions,  respectively,  as  a  function  of  d^  at  Reynolds  numbers  50, 
100,  and  150.  The  total  lift  coefficient,  figure  11  (o),  is  positive  when  the  two  spheres 
are  close  (d,  <  7.9  for  Re  =  50,  d,,  <  4  for  Re  =  100,  and  d,  <  3.4  for  Re  =  150).  That 
is,  the  two  spheres  repel  each  other,  and  the  repulsion  is  stronger  the  closer  they  are. 
Om"  results  show  that  both  the  viscous  and  pressure  contributions  have  an  important 
effect  on  the  repelling  force,  but  the  pressure  contribution  is  more  dominant  when 
Re  ^  100  (compare  figure  ll6,c).  On  the  other  hand,  the  total  lift  coefficient  is 
negative  and  relatively  small  -  that  is,  the  two  spheres  attract  each  other  weakly 
-  at  intermediate  separation  distances  (7.9  <  d^  <  21  for  Re  =  .50,  4  <  d,  <  21  for 
Re  =  100,  and  3.4  <  d,  <  21  for  Re  =  150).  At  these  distances,  the  pressure  is  the 
main  contributor  to  the  attraction  force  at  all  Reynolds  numbers.  The  smaller  the 

FiotTRE  1 1.  Lift  coefficients  of  the  solid  spheres  as  a  function  of  d,  at  Re  =  .50,  100,  and  150 :  (a)  total 
lift  coefficient;  (6)  viscous  contribution  to  lift;  (c)  pressure  contribution  to  lift.  O.  fJe  =  30;  E, 
Re  =  100;  A.  Re  -  150. 
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Figtoe  12.  Moment  coefficient  of  the  solid  spheres  as  a  function  of  d,  at  Se  =  50,  100,  and  150. 

Symbols  as  figure  11. 

Reynolds  number  is,  the  smaller  the  pressure  effect,  the  weaker  the  attraction,  and 
the  narrower  the  region  of  attraction.  When  ^  21,  however,  the  lift  vanishes,  and 
the  two  spheres  have  no  interactions  at  any  Reynolds  numbers. 

As  discussed  in  §3.2.1,  when  the  two  spheres  are  in  close  proximity,  the  front 
stagnation  point  is  shifted  toward  the  x,y  symmetry  plane,  and  the  fluid  particles 
accelerate  faster  in  the  lower  left  region  than  in  the  upper  left  region  of  the  sphere. 
This  difference  in  acceleration  results  in  a  net  negative  pressure  gradient  normal  to 
and  away  from  the  x,y  symmetry  plane,  contributing  to  the  repulsion  between  the- 
two  spheres.  The  shear  stress  is  also  higher  h:  the  lower  left  region  than  in  the  upper 
left  as  shown  in  figure  9.  Furthermore,  owing  to  its  inclination  with  the  x-axis,  the 
shear  force  in  the  lower  left  region  contributes  to  both  the  lift  (parallel  to  the  z-axis) 
and  drag  (parallel  to  the  z-axis),  whereas  the  shear  force  at  the  top  of  the  sphere 
contributes  mainly  to  the  drag.  Therefore,  both  the  pressure  and  shear  forces 
contribute  to  a  positive  lift  force  (i.e.  the  two  spheres  repel  each  other)  when  the  two 
spheres  are  close. 

On  the  other  hand,  when  the  two  spheres  are  in  the  intermediate  separation 
regime,  the  velocity  vector  distributions  show  that  the  front  stagnation  streamline 
is  almost  straight,  and  thus  the  flow  in  the  lower  left  region  is  not  affected  by  the 
presence  of  the  other  sphere.  Nevertheless,  the  gap  between  the  two  spheres  causes 
the  flow  to  accelerate  slightly  faster  on  the  top  of  the  sphere  than  on  the  bottom  and, 
as  a  result,  the  average  pressure  in  the  lower  region  is  slightly  higher  than  that  in  the 
narrow  gap.  Thus,  the  two  spheres  attract  each  other  weakly,  and  the  attraction 
force  is  mainly  due  to  the  pressure  distribution.  The  shear  force,  nearly  parallel  to  the 
z-axis  at  the  top  of  the  sphere,  contributes  mainly  to  the  drag  but  not  to  the  lift. 

3.2.4.  yioment  and  drag  coefficients 

Figure  12  shows  the  moment  coefficient  as  a  function  of  dimensionless  distance  at 
Reynolds  numbers  50,  100,  and  150.  The  moment  coefficient  is  positive  when  the  two 
spheres  are  close  (d,  <  4.6  for  Re  =  50,  d^  <  2.5  for  Re  =  100.  and  d,  <  1.96  for  Re  = 
150),  that  is,  the  two  spheres  experience  positive  torque,  and  the  torque  becomes 
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Fiottre  13.  Drag  coefi&cieat  of  the  solid  spheres  as  a  function  of  at  Re  =  30,  100,  and  150. 

Symbols  as  figure  11. 

stronger  the  closer  they  are.  On  the  other  hand,  the  moment  coefficient  is  negative 
at  intermediate  separation  distances.  The  moment  essentially  vanishes,  and  the  two 
spheres  have  negligible  interactions  with  each  other  when  d^'^21.  In  our 
calculations,  the  solid  spheres  were  not  allowed  to  rotate  under  the  influence  of  the 
torque. 

When  the  two  spheres  are  close,  the  shear  stress,  on  average,  is  higher  aroimd  the 
lower  part  of  the  sphere  than  around  the  top,  and  thus  they  experience  positive 
torque.  On  the  other  hand,  when  the  spheres  are  at  the  intermediate  separation 
distances,  sUghtly  higher  velocity  in  the  gap  leads  to  slightly  higher  shear  stress  on 
the  top  than  on  the  bottom  of  the  sphere,  and  this  causes  the  spheres  to  experience 
weak  negative  torque. 

We  note  that  the  torque  acting  on  the  sphere  is  relatively  small,  and  the  moment 
coefficient  is  less  than  1  %  of  the  drag  coefficient  for  aU  the  separation  distances  and 
Reynolds  numbers.  The  main  reason  for  this  is  that  the  torque  depends  only  on  the 
distribution  of  the  shear  stresses  {r^^  and  and,  as  shown  in  figure  9,  the  shear 
stress  on  the  top  of  the  sphere  counteracts  that  on  the  left  bottom  of  the  sphere. 

Figure  13  shows  the  drag  coefficient  as  a  function  of  the  dimensionless  distance  at 
Reynolds  numbers  50,  100,  and  150.  The  drag  increases  with  decreasing  d^  when  d^ 
is  less  than  4  for  all  Reynolds  numbers.  It  increases  slightly  with  increasing  d^  at 
intermediate  separation  distances,  and  eventually  tends  to  that  of  a  single  sphere 
when  ^  21.  The  drag  increases  as  the  two  spheres  get  close  because  the  shear  stress 
on  the  sphere  is  increased  and  the  pressure  distribution  is  changed  owing  to  the  flow 
acceleration  on  the  lower  left  region  as  well  as  in  the  gap  between  them,  as  shown  in 
figures  8  and  9. 

3.3.  Interactions  of  two  liquid  spheres 

In  the  analysis  of  the  flow  field  past  two  liquid  spheres,  we  use  a  viscosity  ratio 
(internal  fluid  to  external  fluid)  of  25  and  density  ratio  of  300.  These  values  are 
typical  of  liquid-hydrocarbon  fuel  in  a  high-temperature  high-pressure  surrounding 
gas  generally  encountered  in  gas  turbine  combustors  (Raju  «k  Sirignano  1990). 
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FiGTrsE  14.  Velocity  vector  fields  in  the  principal  plane  for  a  liquid  sphere  at  Re  =  100  for 
d,  =  2;  (a)  ertemal  flow;  (6)  internal  flow. 


Figube  13.  Velocity  vector  fields  in  the  principal  plane  for  a  liquid  sphere  at  Re  =  100  for 
d,  =  1.3:  (a)  external  flow;  (6)  internal  flow. 


As  in  the  solid-sphere  case,  we  examine  the  flow  field  for  two  liquid  spheres  in  the 
(x,  2)-plane  of  symmetry,  the  principal  plane,  where  the  narrowest  path  between  the 
liquid  spheres  is  encotmtered.  Figure  5{a-c)  discussed  in  §3.2.1  can  also  represent 
t3rpical  streamlines  in  the  external  flow  of  liquid  spheres.  However,  there  are 
differences  from  the  solid-sphere  case.  First,  the  angle,  measured  from  6  =  0,  at 
which  separation  occurs  on  the  sphere  surface  is  much  higher  than  that  of  the  solid 
sphere.  Second,  a  closer  examination  of  the  velocity  plot  (figures  14a  and  15o)  in  the 
wake  region  indicates  that  the  separating  streamline,  instead  of  being  nearly  normal 
to  the  sphere  surface,  now  curves  closer  to  the  sphere  siirface,  producing  a  ‘  squashed  ’ 
recirculation  zone.  This  behaviour  was  also  seen  in  the  velocity  vector  field  of 
axisymmetric  flow  in  figure  4(a).  The  length  of  the  recirculating  eddy  is  also  slightly 
smaller  than  that  of  the  solid  sphere. 

Figure  14  (a,  b)  shows  the  velocity  vector  fields  of  the  external  and  internal  flows, 
respectively,  in  the  principal  plane  at  Reynolds  number  100  where  the  two  spheres 
are  separated  by  a  distance  d,  =  2.  A  secondary  eddy  in  the  liquid-sphere  stem  region 
is  evident  in  both  the  upper  and  lower  regions  in  the  principal  plane,  but  the  eddy 
centres  in  both  regions  are  asymmetrical  with  respect  to  the  :  =  0  plane.  Also,  these 
eddies  are  concomitant  with  the  occurrence  of  the  eddies  in  both  regions  in  the 
external  flow.  Figure  15  (a,  b)  shows  the  velocity  vector  fields  of  the  external  and 
internal  flows  in  the  principal  plane  at  Reynolds  number  100  for  the  case  of  =  1.5. 
The  secondary  internal  eddy  in  the  liquid-sphere  stem  region  exists  only  in  the  lower 
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Fioxtbe  16.  Total  lift  coefficient)!  of  the  liquid  spheres  as  a  function  of  at  Re  =  50  (O).  100 

(□),  and  150  (A)- 
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Figtoe  17.  Moment  coefficient  of  the  liquid  spheres  as  a  function  of  d,  at  Re  =  50,  100,  and 

150.  Symbols  as  figurs  16. 


region,  and  the  secondary  eddy  in  the  upper  region  no  longer  exists.  The  vanishing 
secondary  internal  eddy  in  the  upper  region  is  accompanied  by  the  disappearance  of 
the  recirculating  eddy  in  the  upper  region  in  the  external  flow. 

Calculations  of  the  lift,  moment,  and  drag  coeflBcients  were  performed  for 
dimensionless  distances  from  the  liquid  sphere  centre  to  the  symmetry  plane  between 
two  liquid  spheres  in  the  range  1.5  ^  d,  ^  25,  for  a  viscosity  ratio  (liquid  to  gas)  of 
25  and  density  ratio  of  300  at  Reynolds  numbers  50,  100,  and  150.  Figures  16,  17,  and 
18  show  the  coefficients  of  total  lift,  moment,  and  drag  as  a  function  of  dimensionless 
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FiGmE  18.  Drag  coefficient  of  the  liquid  spheres  aa  a  function  of  i,  at  Rt  =  50,  100,  and  150. 

Symbols  aa  figure  16. 


distance  at  Reynolds  numbers  50,  100,  and  150.  The  coefficients  of  total  lift,  moment, 
and  drag  are  slightly  smaller  in  absolute  value  than  those  for  the  solid  spheres  at  both 
the  repelling  and  attraction  separation  distances  and  at  all  Reynolds  numbers.  The 
lower  coefficients  of  the  liquid  sphere  are  attributed  to  the  surface  motion  of  the 
liquid  sphere  which  reduces  the  velocity  gradient  and  friction  force.  A  smaller  drag 
coefficient  for  the  liquid  sphere  in  axisymmetric  flow  has  been  also  found  in  earlier 
calculations  (Clift  et  al.  1978). 


4.  Conclusions 

Three-dimensional  flow  interactions  between  two  identical  (solid  or  liquid)  spheres 
which  are  held  fixed,  with  the  line  connecting  the  sphere  centres  normal  to  a  uniform 
stream,  have  been  investigated  at  Reynolds  numbers  50,  100,  and  150  as  a  first  step 
towards  understanding  the  three-dimensional  interactions  with  a  large  concentration 
of  particles. 

First,  the  interactions  between  two  solid  spheres  have  been  investigated  for  a 
dimensionless  distance  in  the  range  1.5  ^  do  ^  25. 

The  two  spheres  repel  each  other  when  they  are  close  (do  <  7.9  for  Re  =  50, 
do  <  4  for  Re  =  100,  and  dg  <  3.4  for  Re  =  150),  and  the  repulsion  is  stronger  the 
closer  thaij  are.  On  the  other  hand,  the  two  spheres  attract  each  other  weakly 
at  intermediate  separation  distances  (7.9  <  d,  <  21  for  Re  =  50,  4  <  do  <  21  for 
Re  =  100,  and  3.4  <  d,  <  21  for  Re  =  150).  For  d,  ^  21,  however,  the  lift  vanishes, 
and  the  two  spheres  do  not  interact  at  any  Reynolds  numbers. 

The  two  spheres  e.xperience  positive  torque  when  they  are  close  (d„  <4.6  for  Re  = 
50,  do  <  2.5  for  i?e  =  100,  and  d,  <  1.96  fori?e  =  150),  and  the  torque  is  stronger  the 
closer  they  are.  On  the  other  hand,  the  moment  coefficient  is  negative  at  intermediate 
separation  distances.  The  moment  vanishes,  and  the  two  spheres  do  not  interact 
when  do  ^  21.  The  drag  on  the  spheres  increases  when  dg  is  less  than  4  for  all 
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Reynolds  numbers.  It  increases  slightly  at  intermediate  separation  distances,  and 
eventually,  tends  to  that  of  a  single  sphere  when  ^  21. 

The  flow  structure  ahead  of  each  sphere  is  such  that  the  streamlines  shift  away 
from  the  X,  y  symmetry  plane  due  to  the  flow  blockage  in  the  gap  between  the  two 
spheres  as  they  approach  the  front  stagnation  region.  Also,  interesting  phenomena 
in  the  near  wake  have  been  discovered  as  the  gap  between  the  two  spheres  decreases. 
When  do  =  2,  the  upper  eddy  is  not  formed  by  the  fluid  separating  on  the  upper 
portion  of  the  sphere,  but  rather  by  the  fluid  turning  around  the  lower  eddy  and 
detached  from  the  sphere.  Furthermore,  when  decreases  to  1.5,  both  the  upper 
eddy  and  downstream  stagnation  point  vanish. 

The  interactions  between  two  liquid  spheres  have  been  also  investigated  for  the 
dimensionless  distance  in  the  range  1.5  ^  dj  ^  25  for  a  viscosity  ratio  of  25  and 
density  ratio  of  300  at  Reynolds  numbers  50,  100,  and  150. 

The  magnitudes  of  the  lift,  torque,  and  drag  on  the  liquid  spheres  are  slightly 
smaller  in  absolute  value  than  those  of  the  solid  spheres  at  ail  the  separation 
distances  and  all  Reynolds  numbers.  The  flow  structure  in  the  external  flow  of  the 
Liquid  spheres  is  quite  similar  to  that  of  the  solid  spheres,  except  that  the  separation 
angle  is  much  higher  than  that  of  the  solid  spheres  and  the  separation  streamline  is 
bent  closer  to  the  sphere  surface  producing  a  ‘squashed’  recirculation  zone. 
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Abstract 

A  numerical  simulation  is  performed  for  the  three- 
dimensional  interaction  of  two  moving  droplets  which 
are  injected  into  an  initially  quiescent  fluid  medium. 
The  pressure  and  velocity  fields  around  each  of  the 
droplet  are  modified  due  to  the  presence  of  the  other 
droplet.  Drag  and  lift  are  therefore  different  from  the 
values  for  em  axisynunetric  flow  around  a  single  iso¬ 
lated  droplet.  The  droplets  decelerate  due  to  the  drag 
Eind  changing  their  direction  of  motion  due  to  the  lift. 
By  choosing  the  origin  of  a  noninertial  reference  frame 
at  the  center  of  mass  of  the  two  droplets,  the  Navier- 
Stokes  equations  are  solved  with  a  noninertial  term  in 
an  iterative  manner.  For  small  initial  separation,  the 
lift  forces  are  repelling,  thereby  increasing  the  separa¬ 
tion.  For  larger  initial  separation,  a  slight  attraction 
occurs. 

N  omenclature 

a'g  dimensional  droplet  radius 

do  initial  ratio  of  distance  between  droplet 

centers  to  droplet  diameter 
dt  instcintaneous  dimensionless  distance 

between  droplet  centers 

It  instantaneous  dimensionless  position 

in  y  direction 

,  N2,  A3  numbers  of  grids  in  y,  (  directions 
Re  initial  Reynolds  number  based  on  droplet 

diameter, 

Rct  instzintaneous  droplet  Reynolds  number, 

Ui^t’laolv 

C/j  ,  initial  dimensional  droplet  velocity 

f/j ,  instantaneous  dimensional  droplet  velocity, 

f  instantaneous  dimensional  droplet  velocity 
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in  x-direction 

lA, ,  instantaneous  dimensional  droplet  velocity 

in  y-direction 

t  time  normailized  by  a'^/U^  „ 

x,y,2  Cartesian  coordinates 

Greek  symbols 

^,T)X  nonorthogonal  generalized  coordinates 

T  viscous  stress  tensor 

Superscript 

'  dimensional  quantity 

Subscript 

o  initial  quantity 

1.  Introduction 

Typical  dbpersed  two-phase  flows  encountered  in  en¬ 
gineering  applications  such  as  combustion  systems  and 
chemical  processes  have  regions  of  large  concentration 
of  particles  or  droplets.  In  such  regions,  the  effect  of 
neighboring  droplets  modifies  the  ambient  conditions  m 
the  flow  near  any  given  droplet.  The  droplet  trajectory 
and  the  coefficients  of  drag,  lift,  and  moment  can  be 
significantly  affected  by  the  modified  flow  field  due  to 
neighboring  droplets. 

The  geometricjd  configurations  of  droplets  in  a  real 
region  of  large  concentration  are  complex  and  subject 
to  uncertainty.  Droplet  arrays  as  discussed  by  Sirig- 
nano  [1],  although  artificial,  can  provide  information 
on  droplets  interaction  and  give  a  detailed  analysis  of 
the  problem. 

Several  investigators  [2-7]  have  performed  research  on 
axisymmetric  configurations  of  interacting  particles  (or 
droplets)  in  the  wake  of  another  particles  (or  droplets). 

Recently,  Kim  ei  al.  [8]  investigated  three- 
dimensional  flow  interactions  with  two  identical  spheri¬ 
cal  droplets  which  were  held  fixed  relative  to  each  other 
in  the  transverse  direction  against  the  uniform  stream 
at  Reynolds  number  0(100);  they  also  studied  two  in- 
dentical  solid  spheres  in  the  same  situation.  They  de¬ 
termined  the  efiects  of  three-dimensional  interactions 
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on  the  lift,  moment,  and  drag  coefficients  as  a  function 
of  the  dimensionless  distance  between  the  two  spheres 
and  Reynolds  number  and  also  discovered  interesting 
near-wake  flow  patterns  as  the  gap  between  the  two 
spheres  decreased. 

In  the  present  paper,  we  extend  the  work  of  Kim  et 
al.  [8]  and  study  the  droplets  interaction  in  a  more  real¬ 
istic  situation  where  two  identical  droplets  are  injected 
and  then  move  side-by-side  into  initially  quiescent  fluid 
medium.  The  droplets  will  be  decelerating  due  to  the 
drag  and  changing  their  direction  of  motion  due  to  the 
lift.  By  placing  the  origin  of  a  noninertial  reference 
frame  at  the  center  of  mass  of  the  two  droplets,  the 
Navier-Stokes  equations  to  be  solved  include  a  nonin¬ 
ertial  term,  which  will  be  evaluated  from  Newton’s  sec¬ 
ond  law  for  the  droplet  motion. 

2.  Formulation  and  numerical  solution 

We  consider  an  unsteady,  three-dimensional,  incom¬ 
pressible,  laminar  flow  generated  by  two  indenticaJ 
(spherical)  droplets  injected  into  initially  quiescent  con¬ 
stant  property  Newtonian  fluid  and  then  moving  side- 
by-side  lying  in  the  same  plane  as  shown  in  Figure  1, 
where  dt  denotes  the  instantaneous  distance,  normal¬ 
ized  by  the  droplet  radius,  from  the  droplet  center  to 
the  y-z  symmetry  plane  between  the  two  droplets.  We 
neglect  the  net  gravity  force  acting  on  the  droplet  and 
also  assume  that  the  Weber  number  is  small  enough 
that  the  droplet  remains  of  spherical  shape.  We  choose 
the  origin  of  a  nonrotating  noninertial  reference  frame 
at  the  center  of  mass  of  the  two  droplets  and  the  angle 
of  initial  droplet  motion  =  0.  Far  upstream,  the 
flow  is  uniform  and  has  an  instantaneous  velocity 
j  parallel  to  the  y-axis.  It  is  noted  that  there  are  two 
symmetry  planes  in  Figure  1  One  is  the  x-y  plane  in 
which  the  centers  of  the  two  droplets  lie,  and  the  other 
is  the  y-z  plane  which  is  midway  between  the  droplet 
centers. 

As  shown  in  Figure  1,  we  utilize  three  coordinate 
systems  in  our  formulation,  the  Cartesian  coordinates 
(x,y,z)  for  the  gas  phase  whose  origin  is  at  the  center 
of  mass  of  the  two  droplets,  the  Cartesian  coordinates 
{xi,yi,zi)  for  the  liquid  phase  whose  origin  is  at  the 
droplet  center,  euid  also  the  nonorthogonal  generalized 
coordinates  (^,?7,C)-  ^  is  the  radial,  r]  is  the  angular, 
and  C  is  the  azimuthal  direction.  Due  to  symmetry 
of  the  geometry,  the  physical  domain  is  chosen  as  one 
quarter  of  an  ellipsoid-like  space. 

For  the  gas  phase  outside  the  droplet,  the  flow  do¬ 
main  is  bounded  by  1  <  ^  <  A'l,  1  <  ;;  <  N2, 
I  <  C  <  N3  +  1.  ^=1  and  correspond  respec¬ 
tively  to  the  droplet  surface  and  the  farfieid  boundary 
surrounding  the  droplet,  rj  =  I  and  Nn  denote  respec¬ 
tively  the  positive  z/-axis  and  the  x;  —  yi  plane  outside 
the  droplet.  C  =  1  ‘"^3  +  1  are  the  same  plane  and 

refer  to  the  xj  —  z;  plane  outside  the  droplet. 

For  the  liquid  phase  inside  the  droplet,  the  domain 


is  1  <  $,  <  Nil,  I  <  11  <  N2,  I  <  Cl  <  N3  +  1.  =  I 

and  Nu  correspond  to  the  center  and  the  surface  of  the 
droplet,  respectively.  r;j  =  1  and  N2  denote  the  positive 
xj-axis  and  the  i/  —  yi  plane,  respectively.  0  =  1  and 
N3  +  I  are  the  same  plane  and  refer  to  the  x/  —  zj 
plcine.  Within  the  droplet,  0  =  constant  are  a  family 
of  concentric  spherical  surfaces.  Uniform  spacing  {6C  = 
St)  =  SC  =  1)  is  used  for  the  generalized  coordinates  in 
both  phases  for  convenience. 

The  nonorthogonal  generalized  coordinate  system 
for  the  external  flow  can  be  easily  adapted  to  three- 
dimensional  arbitrary  geometries.  We  solve  the  time- 
dependent  equations  of  continuity  and  Navier-Stokes 
with  cdl  associated  conditions,  as  will  be  discussed  in 
detail  in  section  2.2. 

2.1  Governing  equations  and  boundary  condi¬ 
tions 

Since  our  goal  is  to  study  the  flow  interaction  with 
the  two  droplets,  we  present  the  equations  for  the  gas 
and  liquid  phases.  The  continuity  and  momentum 
equations  in  both  phases,  and  the  boundary  conditions, 
are  nondimensionalized  using  the  droplet  radius  a'  as 
the  characteristic  length  and  Uj  ^  as  the  characteristic 
velocity. 

Gas  phase 


v*v  =  o 


(1) 


dV, 

dt 


-f_  =  -Vp  + 


(2) 


where  V,  is  the  velocity  of  the  center  of  mass  of  the 
two  droplets. 

Liquid  phase 


V-V,  =  0 


(3) 


dWi 

dt 


-f 


DYi 

Dt 


-  -Vp;  + 


Re. 


-V=V,. 


(4) 


where  is  the  droplet  velocity.  In  these  first  calcu¬ 
lations  presented  herein,  the  inertial  correction  for  the 
liquid  phase  has  been  neglected. 

The  governing  equations  are  transformed  to  the  gen¬ 
eralized  coordinates  (^,p,C)i  which  allow  for  any  three- 
dimensional  body  of  arbitrary  shape.  The  numerical 
integration  of  the  equations  is  performed  using  a  com¬ 
putational  cubic  mesh  with  equal  spacing  ((5,f  =  Si}  - 
SC  =  1). 

Gas/liquid  interface  conditions 

The  conditions  at  the  interface,  ^  =  1  or  .f;  =  .V;,,  are 
based  on  the  the  principle  of  continuity  of  shear  stresses 
and  tangential  velocities.  Since  the  interface  in  our  flow 
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is  always  spherical  (under  the  assumption  of  small  We¬ 
ber  number),  it  is  more  convenient  to  cast  these  con¬ 
ditions  in  terms  of  spherical  coordinates  (r,0,  <p)  whose 


origin  is  at  the  center  of  the  droplet. 

'^r9,a  —  ^,r9,s  i  (h) 

'^r<p,s  —  '^l,rip,a  )  (^) 

^4.,  =  Vij,,  ,  (7) 

V^.a  =  V,,^,a  .  (8) 


dpi  _  dui  _  dvi 


=  0,  wi  =  0 


at  =  1.  (16) 


dpi  _  dui  _  dvi  _  dwi 
dt]i  ~  drji  drji  ~  dr]i 


at  r)i  =  1.  (17) 


dpi  _  dui  _  dvi 
drji  dr]i  diji 


a.t  Tji  =  iVo.  (18) 


where  Trg^a  and  are  respectively  the  shear  stresses 
on  a  positive  r-plzme  in  the  positive  0  and  <j>  direction, 
and  the  subscript  s  denotes  the  surface  of  the  droplet. 
Because  no  fluid  can  cross  the  surface  of  the  droplet, 
the  normal  velocities  at  the  interface  relative  to  surface 
are  zero  in  both  phases.  The  interface  condition  for 
pressure  is  obtedned  from  the  momentum  equation. 

Due  to  the  different  Cartesian  coordinates  used  in 
each  phase,  the  relation  between  the  velocities  on  the 
interface  is  expressed  as  follows. 


u 

—  Ill  +  Uao.t  I 

(9) 

t; 

=  VI  , 

(10) 

w 

=  Wl  , 

(11) 

The  periodic  condition  is  applied  for  ui,vi,wi,  and  pi 
at  the  two  planes  0  =  1  aJid  Wa  -I- 1. 

The  drag,  lift,  and  moment  coefficients  are  evaluated 
in  dimensional  form  as  follows. 


Fp  —  J  -p'rx-eD  dS'  -1-  J  n-r'-eo  dS' 

Fl  =  J  -p'ri'GL  dS'  +  j  n*T''ei  dS' 
M'  =  J  r'  X  r'  dS', 


(19) 


(20) 

(21) 


where  C/oo,«  is  the  droplet  velocity  in  x-direction  due  to 
the  lift  force  acting  on  the  droplet. 

Gas-phase  boundary  conditions 

In  the  following  equations,  u,  v,  and  w  are  the  gas 
velocities  in  the  x,  y,  and  z  direction,  repectively.  p  is 
the  pressure,  amd  the  subscript  I  denotes  liquid  phase. 

p  =  0,  u  =  0,  u  =  Voo  ,,  w  =  Q  at  ^  =  Wi  (12) 

e.xcept  1  =  0. 


where  ep  =  (— stna<i+cosQrJ),  ei  =  (cosorti+sinaj), 
S'  denotes  the  surface  of  the  droplet,  n  is  the  outward 
unit  normal  vector  at  the  surface,  r'  is  the  position 
vector  from  the  center  of  the  droplet,  and  t'  is  the 
viscous  stress  tensor. 

The  repelling  force  is  assumed  positive.  Due  to  the 
geometrical  symmetry,  only  the  z-component  of  mo¬ 
ment  exists,  and  counter-clockwise  direction  is  assumed 
positive. 

The  nondimensional  coefficients  of  drag,  lift,  and  mo¬ 
ment  are  defined  respectively  as 


dp  _  dv 
dx  dx 


dw 

dx 


=  0,  u  =  0 


at  1  =  0.  (13) 


Cd  = 


F' 

Sd 


d,t 


dp  dll  dv  _  dw 
dij  drj  dr]  drj 


at  r/  =  1.  (14) 


Cl  = 


n 


1  '*5 


(23) 


dr] 


dudv 

—  =  -^  =  0,  10  =  0 
dr]  dr] 


at  =  Nn.  (15) 


The  periodic  condition  is  applied  for  u,  v,  w,  and  p  at 
the  two  pleuies  C  =  1  N3  +  1. 


Liquid-phase  boundary  conditions 


Cm 


M'-k 


kp'u'Z 


7ra„ 


(24) 


2.2  Numerical  solution 

We  solve  numerically  the  complete  set  of  Navier- 
Slokes  tind  continuity  equations  in  each  phase,  subject 
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to  the  boundary  conditions  discussed  in  the  previous 
section.  We  have  developed  a  three-dimensional  im¬ 
plicit  finite-difference  algorithm  to  solve  simultaneously 
the  set  of  the  discretized  partied  differential  equations. 

The  governing  equations  are  represented  in  general¬ 
ized  coordinates  and  the  control  volume  for¬ 

mulation  is  used  to  develop  the  finite-difference  equa¬ 
tions  from  them.  The  method  of  solution  employs 
an  Alternating-Direction-Predictor-Corrector  (ADPC) 
scheme  to  solve  the  time-dependent  Navier-Stokes 
equations.  A  pressure  correction  equation  is  employed 
to  satisfy  indirectly  the  continuity  equation. 

The  overall  solution  procedure  is  based  on  a  cyclic 
series  of  estimate-eind-correct  operations.  The  velocity 
components  are  first  calculated  from  the  momentum 
equations  using  the  ADPC  method,  where  the  pressure 
field  at  the  previous  time  step  is  employed.  This  esti¬ 
mate  improves  as  the  overall  iteration  continues.  The 
pressure  correction  is  calculated  from  the  pressure  cor¬ 
rection  equation  using  the  SOR  successive  overrelaixtion 
method.  The  new  estimates  of  pressure  and  velocities 
are  then  obtained. 

The  szime  procedure  is  performed  in  the  liquid  phase. 
The  governing  equations  of  motion  in  each  phase  are 
solved  in  an  interactive  sequence  through  the  interface 
boundary  conditions  until  convergence  is  achieved  for 
each  time  step  of  the  calculation.  Changes  in  droplet 
velocity  are  also  determined  by  resolving  lift  and  drag 
forces  on  the  droplet  and  applying  Newton’s  second 
law.  This  process  continues  until  the  solution  converges 
at  each  time  step. 

In  the  overall  procedure,  the  sequential  solutions  of 
governing  equations  and  boundary  conditions  with  grid 
and  relative  velocity  adjustment  are  iterated  until  con¬ 
vergence  is  achieved.  After  convergence  is  reached,  the 
drag,  lift,  and  moment  coefficients  are  evaluated  at  the 
prescribed  time  interved. 

Now,  the  moving  droplets  change  their  directions  due 
to  the  interactions  with  the  surrounding  fluid.  How¬ 
ever,  in  order  to  insure  its  accuracy,  the  computational 
grid  must  satisfy  the  condition  that  the  lift  and  torque 
are  zero  for  a  single  droplet  moving  in  any  direction. 
Figures  2(a)  and  2(b)  show  two  different  grid  distribu¬ 
tions  on  the  droplet  surface.  The  grid  in  Figure  2(a)  for 
a  single  droplet  gives  at  finite  angle  Og  nonzero  lift  and 
torque  which  have  the  same  order  of  magnitude  as  those 
of  droplets  separated  by  dg  <  2.  On  the  other  hand, 
the  grid  in  Figure  2(b)  gives  zero  lift  and  torque  at  any 
angle.  We  generate  this  three-dimensional  grid  alge¬ 
braically  by  using  a  set  of  ellipsoids,  cones,  and  planes, 
where  grid  density  is  controlled  by  the  stretching  func¬ 
tion  developed  by  Vinokur  [9].  Figure  3  presents  the 
grid  on  the  x-y  symmetry  plane  with  the  droplet  cen¬ 
ters  at  dimensionless  distance  dg  =  9.  Our  numeri¬ 
cal  solution  scheme  has  been  developed  for  arbitrary 
grid  aystems  so  that  it  is  necessary  only  to  change  the 
subroutines  of  grid  generation  for  flow  geometries  that 


differ  from  ours. 

3.  Results  and  discussion 

We  first  tested  the  three-dimensional  code  by  solving 
the  steady  axisymmetric  flow  past  a  single  droplet  at 
Reynolds  number  100.  Since  the  code  developed  solves 
for  three  Cartesicui  components  of  velocity  in  the  trans¬ 
formed  grid,  an  axisymmetric  test  calculation  still  ex¬ 
ercises  the  fully  three-dimensionzd  aspects  of  the  code. 

3.1  Flow  over  a  single  droplet 

Here  we  discuss  the  flow  generated  by  an  impulsively 
started  single  droplet  into  an  initially  quiescent  fluid 
using  the  three-dimensional  solution  procedure  for  vis¬ 
cosity  ratio  of  25  and  a  density  ratio  of  300  (liquid  to 
gas)  at  Reynolds  number  100.  The  time-dependent  so¬ 
lution  converges  asymptotically  to  the  steady-state  re¬ 
sults.  Table  1  lists  the  drag  coefficients  (Cdp  and  Cdv 
are  respectively  the  pressure  and  viscous  parts  of  C/j) 
as  a  function  of  the  computational  grid  density  which 
are  in  good  agreement  with  the  correlation  of  Rivkind 
&  Ryskin  [10].  Table  2  shows  two  angles  measured 
from  the  front  stagnation  point  where  is  the  angle 
at  which  the  surface  vorticity  changes  its  sign  and  0^  is 
the  angle  at  which  the  surface  velocity  is  zero. 

The  calculations  were  performed  for  three  different 
grids,  (Ni  X  Ni  x  N3  and  Nu  x  N21  x  A^a;)  =  (20  x  11  x  32 
and  10  X  11  X  32),  (30  x  15  x  48  and  15  x  15  x  48  ), 
and  (40  x  21  x  64  and  20  x  21  x  64),  in  .a  computational 
domain  having  an  outer  boundary  located  at  21  droplet 
radii  from  the  droplet  center.  The  solutions  from  the 
three  different  grids  were  stable  and  smooth,  and  each 
takes  dimensionless  times  of  20  to  reach  steady  state. 

Tables  1  and  2  show  that  the  results  of  the  30  x  15  x  48 
grid  differ  by  only  0.8  %  in  the  drag  coefficient  and  0.7 
%  in  the  sepeu’ation  angle  from  those  of  the  40  x  21  x  64 
grid.  Thus,  for  computationzJ  economy  (a  long  time 
period  is  needed  for  our  time-dependent  solution  of  the 
moving  droplets),  we  selected  the  medium  size  grid  30  x 
15  X  48,  and  15  x  15  x  48  in  the  liquid  phase  for  the 
remaining  calculations. 

3.2  Interaction  of  two  moving  droplets 

Calculations  were  performed  for  initial  dimension¬ 
less  distance  dg  =  2  and  9  at  initial  Reynolds  num¬ 
ber  100  for  the  time  period  0  <  t  <  250.  This  time 
period  is  about  one  quarter  of  the  lifetime  for  a  vapor¬ 
izing  droplet  injected  into  a  combustor,  and  so  major 
changes  in  transport  processes  and  flow  field  would  be 
expected  to  occur.  The  viscosity  and  density  ratios  (liq¬ 
uid  to  gas)  for  the  droplet  is  25  and  300,  respectively, 
which  are  typical  of  liquid-hydrocarbon  fuel  in  high- 
temperature,  high-pressure  surrounding  gas  generally 
encountered  in  gas  turbine  combustors. 

Figure  4(a)  shows  the  trajectory  of  one  droplet  for 
dg  =  2,  where  the  numbers  on  the  lines  denote  the 
dimensionless  time.  The  final  location  of  the  droplet 
at  1  =  250  is  (dt,/,)  =  (2.505,  212.6).  Now,  since  the 
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final  location  of  a  single  droplet  at  t  =  250  in  the  case 
of  no  drag  and  lift  forces  acting  on  it  would  be  (dtjt) 
=  (2,  250),  the  figure  indicates  that  the  two  droplets 
are  repelling  each  other  as  well  as  decelerating,  and 
the  change  of  distance  (normalized  by  droplet  radius) 
is  much  higher  (37.4  vs.  0.505)  due  to  the  drag  than 
due  to  the  lift  force.  Figure  4(b)  shows  the  trajectory  of 
one  droplet  for  do  =  9.  The  final  location  of  the  droplet 
at  t  =  250  is  =  (8.968,  213.1).  The  theoretical 

final  location  of  a  single  droplet  at  t  =  250  in  the  case 
of  no  drag  and  lift  forces  would  be  {dtJt)  =  (9>  250); 
therefore,  Figure  4(b)  identifies  that  the  two  droplets 
are  weakly  attracting  each  other  as  well  as  decelerating, 
cind  the  change  of  distance  is  much  higher  due  to  the 
drag  than  due  to  the  lift  force. 

Figure  5  shows  the  lift  coefficients  of  the  droplet  for 
dg  =  2  and  9  as  a  function  of  instantaneous  Reynolds 
number,  where  the  repelling  force  is  taken  as  positive. 
The  lift  coefficient  for  do  =  2  is  positive  during  the  time 
period  0  <  1  <  250,  and  becomes  gradually  smaller  in 
time  because  the  distance  between  the  droplets  is  in¬ 
creasing  2ind  their  directions  of  motion  are  changing  due 
to  the  repelling  force.  Qn  the  other  hand,  the  lift  coeffi¬ 
cient  for  do  =  9  is  negative  during  that  time  period  but 
slowly  goes  towards  zero.  The  change  of  the  distance 
between  the  two  droplets  for  dg  =  9  case  is  negligi¬ 
bly  small  as  shown  in  Figure  4(b).  Thus,  this  result 
also  indicates  that  the  lift  coefficient  slowly  approaches 
zero  when  Reynolds  number  decreases  with  the  fixed 
dimensionless  distance  9.  The  figure  also  shows  that 
rapid  change  occurs  initially  due  to  impulsive  start  of 
the  droplets  and  the  quasi-steady  state  occurs  when 
t  >  30. 

Figure  6  shows  the  moment  coefficients  of  the  droplet 
for  dg  =  2  and  9  as  a  function  of  instcintaneous 
Reynolds  number,  where  counter-clockwise  direction  is 
take  as  positive.  The  moment  coefficient  for  d*  =  2  is 
positive  initially  and  gets  gradually  smaller  and  changes 
its  sign  from  positive  to  negative  at  t  =  150.  On  the 
other  hand,  the  moment  coefficient  for  d,  =  9  is  neg¬ 
ative  during  the  time  period  0  <  f  <  250  but  slowly 
approaches  zero. 

Figure  7  shows  the  drag  coefficients  of  the  droplet  for 
dg  =  2  and  9  as  a  function  of  instantaneous  Reynolds 
number.  The  drag  coefficient  for  d^  =  2  is  higher  than 
that  for  dg  =  9.  In  earlier  time,  the  difference  is  greater 
but  becomes  gradually  smaller  as  time  goes  on.  The 
reason  is  that  the  droplets  for  dg  —  2  are  repelling  each 
other  and  the  distance  between  them  is  increasing,  and 
so  the  drag  becomes  gradually  smaller.  The  drag  co¬ 
efficient  for  do  =  9  is  almost  identiczd  to  that  for  a 
single  droplet  (slightly  higher  that  for  a  single  dronlet 
when  t  >  30.)  Tsuji  et  al.  [10]  amd  Kim  ei  al  uso 
found  that  for  the  fi.\ed  spheres  and  dg  >  4,  the  drag 
coefficient  is  almost  identical  to  that  of  a  single  sphere. 
Figure  7  also  shows  large  drag  coefficient  initially  which 
indicates  that  large  shear  stress  and  pressure  occur  ini¬ 


tially  due  to  impulsive  start. 

Figure  8  shows  the  speed  of  the  droplet,  as  a  func¬ 
tion  of  dimensionless  time.  There  is  a  sudden  drop 
in  the  speed  at  the  moment  of  droplet  injection  and 
then  monotonically  decreases  with  time.  The  speed  for 
do  =  2  is  lower  than  dg  =  9  due  to  the  larger  drag  of 
do  =2. 

Figures  9(a)  and  9(b)  show  the  velocity  vector  fields 
of  gas  and  liquid  phases  at  Ret  =  96.86,  dt  =  2.002, 
and  t  =  19.9  in  the  x-y  plane  of  symmetry,  where  the 
narrowest  path  between  the  two  droplets  is  encoun¬ 
tered.  The  free  stream  is  coming  from  the  left  in  the 
figures.  The  velocity  vectors  were  first  calculated  in 
the  Cartesian  coordinates  fixed  at  the  center  of  mass  of 
the  two  droplets  and  then  converted  in  the  Cartesian 
coordinates  fixed  at  the  center  of  the  droplet.  Due  to 
the  divergence  of  the  front  dividing  streamline  and  the 
shift  of  the  front  stagnation  point  toward  the  y-z  sym¬ 
metry  plane,  the  fluid  particles  accelerate  more  in  the 
lower  left  region  around  the  droplet  than  in  the  upper 
left.  Furthermore,  the  pressure  drop  in  the  lower  left 
region  is  higher  than  the  pressure  drop  through  the  nar¬ 
row  path  between  the  two  droplets.  Shear  stresses  on 
the  lower  left  region  are  also  higher  and  contribute  to 
lift  as  well  as  drag  due  to  their  incUnation  with  y-axis, 
whereas  the  shear  stresses  at  the  top  of  the  droplet  are 
higher  due  to  narrow  path  but  contribute  mainly  to  the 
drag.  As  a  consequence,  the  two  droplets  ait  repelling 
each  other.  These  phenomena  are  explained  in  detail 
in  the  study  of  Kim  et  al.  [8]  about  the  flow  past  fixed 
spheres. 

Figures  10(a)  and  10(b)  show  the  velocity  vector 
fields  of  gas  and  liquid  phases  at  Ret  =  96.94,  dt  =  Sj 
and  t  =  19.9  in  the  x-y  plane  of  symmetry.  The  veloc¬ 
ity  vectors  near  the  front  stagnation  point  indicate  that 
the  front  dividing  streamline  is  similar  to  that  of  a  sin¬ 
gle  droplet  and  the  front  stagnation  point  moves  very 
slightly  away  from  the  y-z  symmetry  plane.  Slightly 
lower  pressure  occurred  in  the  upper  region  than  in  the 
lower  region  due  to  the  relatively  narrow  path  between 
the  two  droplets.  As  a  consequence,  the  two  droplets 
are  weakly  attracting  each  other. 

4.  Conclusions 

A  numerical  simulation  wcis  performed  for  the  three- 
dimensional  interaction  of  two  moving  droplets  which 
are  injected  into  an  initially  quiescent  fluid  medium. 
Calculations  were  performed  for  initial  dimensionless 
distance  dg  —  2  and  9  at  initial  Reynolds  number  100 
for  the  time  period  0  <  1  <  250. 

The  droplets  for  dg  =  2  repel  each  other  as  well  as 
decelerate.  Their  lift  coefficient  is  positive  during  the 
time  period  0  <  I  <  250  and  gets  gradually  smaller  due 
to  the  change  in  the  distance  and  direction  of  motion. 

The  droplets  for  d,  =  9  weakly  attract  each  other 
as  well  aa  decelerate.  Their  lift  coefficient  is  negative 
during  the  time  period  0  <  1  <  250,  but  goes  slowly 
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towards  zero.  This  indicates  that  the  droplets  for  that 
separation  distance  are  tending  towcirds  weaker  repul¬ 
sion  as  the  Reynolds  number  decreases. 

The  drag  coefficient  for  do  =  2  is  higher  than  that 
for  do  =  9.  In  earlier  time,  the  difference  is  greater  but 
becomes  gradually  smaller  as  time  goes  on.  The  drag 
coefficient  for  do  =  9  is  almost  identical  to  that  for  a 
single  droplet. 
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Gas  Liquid 

iVj  X  iVj  X  iV’3  M\i  X  N'it  X  Nzt  Cop  Cqv  Cq  C’q 


T 

i' 


20  X  11  X  32 

10  X  11  X  32 

0.542 

0.586 

1.128 

30  X  13  X  48 

15  X  15  X  48 

0.520 

0.573 

1.092 

40  X  21  X  64 

20  X  21  X  64 

0.511 

0.571 

1.081  1.08 

Table  1.  Drag  coefficients  as  a  function  of  grid  density  at  Re  =  ICO, 
where  *  denotes  the  data  from  the  correlation, 
of  Rivkind  k  Ryskin  (1976) 


Gas 

Liquid 

iVj  X  fV2  X  Nz 

Nu  X  Nztx  Nzi 

^IJ 

9^ 

20  X  11  X  32 

10  X  11  X  32 

126.31 

146.50 

30  X  15  X  48 

15  X  15  X  48 

127.70 

150.13 

40  X  21  X  64 

20  X  21  X  64 

127.74 

151.20 

Table  2.  An^es  at  which  the  surface  vorticity  changes  its  sign  and 
the  surface  velocity  is  zero  as  a  function  of  grid  density 
at  Re  =  100,  where  the  angles  are  measured  from 
the  front  stagnation  point. 


Figure  1.  Flow  geometry  and  coordinates. 


,y 


Figure  2.  Two  different  grid  distributions  on  the  droplet  surface. 


Figure  3.  Grtd  disthbutioa  at  dimeosioolesa  distance  9. 


Figure  6.  Moment  coefficients  as  a  function  of  instantaneous  Rcvnolds 
for  do  =  2  and  9,  where  the  numbers  on  the  curvw  denote 
the  dimensionless  time. 


Figure  I .  Drag  coefficrents  as  a  funccioa  of  lostantaiieous  Reynolds  number 
<^0  =  2  and  9,  where  the  numbers  on  the  curves  denote 
the  dimensionless  time. 
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Figure  9(a)  Velocity  vector  field  of  the  gas  phase  in  x-y  plane 
at  Re,  =  96.36,  d,  =  2.002,  and  (  =  19.9 


Figure  3.  Speeds  of  the  droplet  as  a  function  of  dimensionless  time 
for  d.  =  2  and  9. 
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Figure  9(b)  Velodcy  vector  field  of  the  liquid  phase  in  x-y  plane 
at  Re,  =  96.36.  d,  =  2.0U2,  and  (  =  ’.9.9 
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Figure  10(a)  Velocity  \-ector  field  of  the  gas  phase  in  x-y  plane 
Re,  =  96.94,  d,  =  9,  and  J  =  19.9 


Figure  10(b)  Velocity  vector  field  of  the  liquid  phase  in  .x-y  plane 
Re,  =  96.94,  d,  =  9,  and  t  =  19.9 
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Abstract 

The  present  study  extends  the  previous  droplet  models  [1,  2]  to  investi¬ 
gate  numerically  the  system  of  three  droplets  vrhich  are  moving  in  tandem 
with  respect  to  the  free  flow.  The  purposes  of  this  study  are  to  study  the 
wake  efi'ect  of  the  lead  droplet  on  the  downstream  droplets  and  to  exam¬ 
ine  the  effects  of  initial  spacing  on  the  total  system.  The  effects  of  variable 
thermophysical  properties,  transient  heating  and  internal  circulation  of  liq¬ 
uid,  deceleration  of  the  flow  due  to  the  drag  of  the  droplet,  boundary-layer 
blowing,  and  moving  interfaice  due  to  surface  regression  as  well  as  relative 
droplet  motion  are  included.  The  results  are  compared  with  those  of  an  iso¬ 
lated  droplet  [1]  as  well  as  those  of  the  two-droplet  system  [2]  to  investigate 
the  effect  of  the  presence  of  the  third  droplet.  The  intereu:tion  effects  from 
the  downstream  or  upstream  droplet  are  identified.  The  transport  rates  of 
droplets  are  reduced  from  the  values  for  an  isolated  droplet,  and  values  for 
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the  downstream  droplets  axe  profoundly  less  than  those  for  the  lead  droplet. 
The  difference  in  transport  rates  is  large  between  the  first  two  droplets;  how¬ 
ever,  it  becomes  insignificant  between  the  second  and  the  third  dropIets.The 
modifications  to  the  transfer  correlations  for  an  isolated  droplet  needed  to 
accoimt  for  the  interaction  effects  are  determined. 


Nomenclature 

Bh  CpgjiirJT'oo  -  effective  heat  transfer  number 

specific  heat  of  gas  phase 

D\2  non-dimensional  droplet  spacing 

between  the  first  and  the  second  droplets 
Z)23  D2Z*/R\  q,  non-dimensional  droplet  spacing 

between  the  second  and  the  third  droplets 
L  Uf{T^Cp\  q),  latent  heat  of  vaporization 

Nu  ^  ^  Nusseit  number 

R  RfRfi  Q,  non-dimensional  instantaneous  droplet  radius 

Q*  heat  flux 

gas-phase  Reynolds  number 

Sh  f*?*^^**”^*^^  Sherwood  number 

s.p.  separation  point  at  droplet  surface  (degree) 

T  T' I temperature 

t'  time 

Uoo  instantaneous  free-stream  velocity 

Yi  mass  fraction 


Greek 


>^g  >^gl'^g,oo^  conductivity  of  gas  phase 

f^e  l^gft^g,oo  viscosity  of  gas  phase 

pfg  density  of  gas  phase 

'’’ifa  S^s  hydrodynamic  diffusion  time 

Superscript 
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'  dimensional  quantity 

Subscripts 
/  fuel 

film  film  conditions,  average  of  ambient  and  surface  conditions 
g  gas  phase 

/  liquid  phase 

*  numerical  index  for  droplets,  l=lead  droplet,  2=second  droplet, 

3=third  droplet 
iso  isolated  status 

s  average  surface  condition 

0  initial  condition 

oo  firee  stream  condition 


1  Introduction 

In  realistic  spray  situations,  the  fuel  is  usually  introduced  into  the  combustor 
as  a  stream  of  liquid  that  breaks  into  droplets.  The  droplets  subsequently 
vaporize  in  the  convective  gas  stream  to  form  the  air-fuel  mixture.  Typi¬ 
cally,  the  fuel  is  of  su£5ciently  low  volatility  that  vaporization  is  an  impor¬ 
tant  controlling  factor  in  the  estimation  of  combustion  rates.  Usuadly,  in  the 
dense-spray  regions  such  as  regions  near  the  fuel  nozzle,  the  droplet  spacing 
is  so  small  that  the  interaction  effects  would  modify  the  droplet  behavior 
significantly.  In  order  to  obtain  the  qualitative  modification  of  the  transfer 
coefficients  in  the  practical  dense  spray  calculation,  it  is  then  necessary  to 
consider  the  interactions  among  a  stream  of  multiple  moving  droplets.  How¬ 
ever,  a  detailed  and  accurate  simulation  of  hundreds  of  droplets  will  be  very 
time  consuming  and  very  difficult  to  perform.  In  fact,  the  behavior  of  trmling 
droplets,  which  follow  the  first  two  or  three  droplets,  can  be  estimated  from 
that  of  the  first  two  or  three  droplets  on  account  of  the  |>eriodical  nature 
of  linear  droplet  arrangements.  This  research  addresses  the  inter<iction  of 
three  vaporizing  droplets  moving  coUinearly  which  represents  a  model  of  an 
injected  stream  of  fuel  droplets. 
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There  is  a  lack  of  a  detailed  investigation  of  multiple-droplet-interactions 
involving  variable  properties,  transient  heating  and  internal  circulation  of 
droplets  in  the  literature.  Even  the  simplified  numerical  computations  of 
three-droplet  dynaunics  are  rarely  found.  TaJ  et  al.  [4]  used  a  multisphere 
cylindrical  cell  model  to  study  hydrodynamics  and  heat  transfer  in  assem¬ 
blages  of  spheres.  They  found  the  hydrodynamic  solution  and  Nusselt  num¬ 
ber,  defined  by  using  average  bulk  temperature  of  the  cell  unit  inlet,  to  be 
perodic.  Tong  and  Chen  [5]  have  included  vaporization  in  Tal’s  model  to  in¬ 
vestigate  the  effects  of  droplet  spacing  on  heat  and  mass  transfer  of  droplets 
in  a  liquid  droplet  array.  A  three-droplet  auray  in  a  cylindrical  duct  has  been 
used  to  obtain  correlations  between  Nusselt  number  and  local  ambient  prop¬ 
erties  for  each  droplet.  The  cylindrical  cell  model  is  somewhat  idealized  and 
some  aissumptions  must  be  imposed  on  the  cell  boundaries.  Hence,  their  re¬ 
sults  must  be  verified  by  calculations  fiom  an  advanced  model.  Kleinstreuer 
et  al.  [6]  used  a  finite-element  microscale  analysis  to  find  the  drag  coefficients 
of  interacting  spheres  in  a  linear  array  and  a  boundary-layer  analysis  for  va¬ 
porizing  droplets  to  simulate  coupled  transfer  processes  for  three  interacting 
droplets  in  a  one-dimensional  trajectory.  Their  solution  is  basically  the  com¬ 
bination  of  three  single-droplet-solutions  with  an  “effective  approach  strekn 
temperature”  to  account  for  interactions.  The  heat  and  momentum  trans¬ 
fer  of  spheres  in  a  linear  array  was  studied  by  Tsai  and  Sterling  [7]  where 
a  finite  difference  method  combined  with  a  body  fitted  grid  was  employed. 
The  transient  effect  and  mass  transfer  was  excluded  in  this  numerical  calcu¬ 
lation.  The  results  of  Nusselt  number  and  drag  coefficients  of  droplets  are 
qualitatively  in  agreement  with  the  findings  of  Tong  and  Chen  [5]  and  Tal  et 
>1.  [4). 

In  the  present  study,  the  momentum,  heat  and  mass  transfer  of  three  in¬ 
teracting,  vaporizing  droplets  are  taken  into  account  and  the  relative  motion 
due  to  different  drag  forces  that  the  droplets  have  experienced  is  included. 
We  aim  to  understand  the  Wcike  effects  on  the  trzmsport  rates  of  the  down¬ 
stream  droplets.  The  primary  emphasis  will  be  placed  upon  the  interaction 
effects  due  to  different  initial  spacings. 

The  schematic  flow  configuration  is  presented  in  Figure  1  where  the 
flow  passing  over  three  vaporizing  droplets  moving  in  tandem  is  shown.  The 
flow  is  laminar  and  aocisymmetric  with  initially  uniform  ambient  conditions 
specified  by  ^/.oo  =  0.  The  initial  droplet  spacing  is 

also  prescribed.  Dij  is  the  non-dimensional  droplet  spacing  (with  respect 


4 


to  the  initial  radius  of  the  first  droplet)  between  the  droplet  and  the 
droplet.  The  frame  of  reference  is  fixed  to  the  center  of  the  lead  droplet  and 
an  jwrcount  will  be  made  for  this  noninertid  frame  by  the  use  of  a  reversed 
pseudo  force.  The  problem  can  be  viewed  as  an  impulsively-started  flow  over 
a  fixed  droplet  and  two  moving  downstream  droplets  aligned  in  tandem. 

The  full  consideration  of  forced  convection  of  the  gas  phase,  transient 
deceleration  of  the  flowfield  due  to  the  drag  force,  internal  circulation  and 
transient  heating  of  the  liquid  phase,  variable  properties,  and  transport  pro¬ 
cesses  occurring  at  the  vaporizing  droplet  interface  required  to  solve  the 
Navier-Stokes  equations,  energy  equation  and  species  equation,  combined 
with  appropriate  boundary  conditions  simultaneously.  Also,  in  order  to  con¬ 
sider  the  moving  boundaries  due  to  surface  regression  and  relative  droplet 
motion,  a  general  method  of  generating  boundany-fitted  coordinate  systems 
is  required. 

The  aodsymmetric  governing-equations  in  cylindricad  coordinates  amd  their 
corresponding  finite-difference  equations,  amd  numerical  procedures  are  given 
in  Chiang  [3].  The  computer  codes  employed  in  our  previous  research  [1,  2] 
have  been  modified  to  dead  with  the  present  three-droplet  am’angement  The 
modification  involves  the  grid  generation  routine  as  well  as  the  routines  which 
handle  the  relative  motion  and  spacings  between  droplets.  The  cadling  se¬ 
quence  amd  pau’ameter  transfer  aunong  subroutines  have  been  auiapted  to  ac¬ 
commodate  new  vairiables  and  the  increase  in  memory  size. 


2  Results  and  Discussion 

Six  production  runs  simulating  three  equi-sized  droplets,  with  the  same  initial 
droplet  Reynolds  number  but  with  different  initiad  droplet  spacings,  moving 
collineau-ly  have  been  performed.  The  values  of  physicad  pamameters  employed 
in  each  case  ame  the  same  as  the  base  case  of  the  two-interaicting-droplet 
study  [2].  The  vadues  of  initiad  spacings  in  each  caise  are  given  in  Table  1. 

The  computations  are  performed  on  a  CRAY  Y-MP  supercomputer.  Fig¬ 
ure  1  adso  shows  the  typical  grid  distributions  at  the  beginning  amd  at  the 
final  computational  time  for  the  case  of  droplet  coalescence.  It  is  noted  that 
the  computation  is  stopped  when  the  droplet  spacing  is  reduced  below  2.6 
or  whenever  the  downstream  droplet  approaches  the  outer  computational 
boundary,  since  the  grid  generation  routine  may  generate  an  overskewed  grid 
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system  under  these  conditions.  Often  the  computation  terminated  at  a  very 
early  stage  of  the  droplet  lifetime  during  which  most  of  heat  flux  to  the 
droplet  surfzice  would  be  transferred  to  the  interior  of  the  droplet.  Hence, 
the  effective  transfer  number  {Bh  =  ^fg^jurJXoo  ~  ^j)(1  “  is  quite 

small.  The  aunount  of  mass  evaporated  is  insignificant,  as  a  result  the  evap- 
oration  rates  of  droplets  are  not  reported  here.  The  time  scale  used  in  the 
following  discussion  is  the  gas-phase  hydrodynamic  diffusion  time  scale. 

The  representative  results  to  characterize  the  droplet  behavior  are  sum¬ 
marized  below. 

The  global  contours  of  results  from  Case  1  (with  D12=D23=12)  and  Case 
4  (D12=D23=6)  are  compared  in  order  to  study  the  effect  of  initial  spacing  on 
the  flow  field.  The  vorticity  distributions  and  isotherms  for  the  gas  and  liquid 
phases  are  presented  in  the  top  and  bottom  portions  of  each  plot  of  Figure  2, 
respectively.  The  results  of  Case  1  at  two  different  times  are  portrayed  in  the 
upper  two  plots.  At  the  early  time,  droplets  experience  the  convective  effect 
such  that  high  vorticity  gradients  occur  at  the  front  portion  of  the  droplet, 
and  an  asymmetric  distribution  is  developed.  The  flow  field  surrounding  each 
droplet  resembles  that  surrounding  an  isolated  droplet  [1].  At  five  gas-phase 
diffusion  times,  the  D12  is  reduced  from  12  to  2.96.  The  vorticity  convected 
downstream  from  the  first  droplet  has  directly  touched  the  second  droplet 
and  shifted  its  vortex  center  to  the  equatorial  plane  since  less  convection 
has  acted  on  the  downstream  droplet.  The  second  droplet  is  well  protected 
by  the  vorticity  wake  of  the  lead  droplet,  while  the  third  droplet,  spaced 
9.5  droplet  radii  away  from  the  second  droplet,  behaves  qu2ditatively  as  an 
isolated  droplet. 

For  the  case  of  small  initid  droplet  spacing  (Case  4),  the  interaction 
effects  are  expected  to  be  strong.  The  velocities  approaching  the  second 
and  the  third  droplets  are  more  or  less  similau'  but  are  considerably  smaller 
than  that  appro2iching  the  lead  droplet.  As  a  result,  the  convective  effect  on 
transport  is  subdued.  In  both  cases,  the  vorticity  distribution  for  the  lead 
droplet  is  distorted  by  the  approaching  of  the  second  droplet. 

The  isotherms  demonstrate  thermal-tramsport  interactions  among  three 
droplets  and  the  gas  phase.  For  the  case  with  sufficient  large  droplet  spacing, 
the  three  droplets  show  the  same  isolated-droplet-isothermal-pattem  at  the 
early  time  as  indicated  in  the  uppermost  plot.  The  third  trailing  droplet 
always  exhibits  the  same  temperature  contour  as  that  of  an  isolated  droplet. 


6 


even  though  the  approax:hing  temperature  and  velocity  have  decreased  from 
the  values  upstream  of  the  lead  droplet.  The  case  with  small  initi2j  spacing 
shows  quite  different  thermal  transport  mechanisms.  The  temperature  gradi¬ 
ents  for  the  downstream  droplets  are  smaller  than  those  for  the  lead  droplet 
due  to  the  action  of  the  thermal  wake.  The  thermal  boundary  layer  thickness 
increases  in  the  downstream  direction  jJong  the  surface  for  the  lead  droplet. 
This  behavior  is  opposite  for  the  third  trailing  droplet  in  the  close  spacing 
case. 

The  transport  properties  of  trailing  droplets  are  substantially  affected  by 
the  fuel  vapor  convected  downstream.  The  composition  of  the  mixture  and 
hence  the  density  of  gas-phase  surrounding  each  droplet  will  be  different  for 
each  droplet  depending  on  the  droplet  spacing.  The  mass-fraction  contour 
of  Case  1  at  two  different  times  are  presented  in  Figure  3.  An  envelope 
of  vapor  wake  created  by  the  leading  droplet  and  covering  the  downstream 
droplet  is  identified  even  though  the  droplet  spacings  are  IcU'ge.  This  envelope 
serves  to  reduce  the  exchanges  of  momentum,  mass  and  energy  between  the 
ambience  and  the  trailing  droplets.  As  time  progresses,  the  droplet  spacing 
from  the  leading  droplet  is  reduced,  the  change  in  spatiad  vapor  distribution 
of  downstream  droplet  is  more  profound. 

The  variations  of  Nusselt  number  along  the  gas/liquid  interface  of  the 
droplets  for  Case  1  and  Case  4  are  shown  in  Figure  4.  In  Case  1,  the  decrease 
of  Nusselt  number  at  the  front  stagnation  portions  of  the  trailing  droplets  are 
caused  by  the  cold  fuel/air  mixture  convected  downstream  from  the  upstream 
droplets.  The  qualitative  variations  of  Nusselt  numbers  of  the  second  droplet 
and  the  third  droplet  in  Case  4  are  totally  different  from  those  in  Case  1. 
The  details  have  been  discussed  in  Chiang(3].  The  third  droplet  shows  the 
same  Nusselt  number  distribution  as  for  the  second  droplet  except  at  the 
separation  region,  where  the  third  droplet  possesses  a  higher  value  of  Nusselt 
number  due  to  the  angular  diffusion  of  heat  flux  from  the  free  stream.  Similzir 
behavior  for  Sherwood  numbers  are  observed  for  the  three  droplets. 

The  surface  shear  stress  distributions  are  shown  in  Figure  5.  The  cascade 
effect  of  upstream  wake  decreases  the  surface  stress  of  downstream  droplets. 
It  is  well  known  that  the  momentum  is  dissipated  along  the  downstream 
direction  due  to  the  wake  action.  Usually  the  recirculating  wake  from  the 
lead  droplet  has  dissipated  a  significant  portion  of  the  momentum  upstream 
of  the  second  droplet.  As  a  result,  the  further  reduction  in  shear  stress  from 
the  second  to  the  third  droplet  becomes  very  limited.  Hence,  the  stress 
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difference  between  the  first  two  droplets  is  much  larger  than  that  between 
two  downstream  droplets  for  Case  1.  Note  that  in  Case  4,  the  shear  stress  for 
the  third  droplet  is  higher  thjin  that  for  the  second  droplet  since  the  second 
droplet  interacts  with  both  neighboring  droplets.  The  negative  shear  stresses 
are  primarily  caused  by  the  contacting  wakes  of  close  neighboring  droplets. 

Figures  6  and  7  present  the  time  variations  of  drag  coefficients  emd  Nus- 
selt  numbers,  respectively.  As  we  expect,  the  transport  rates  of  droplets  are 
reduced  from  the  values  for  an  isolated  droplet,  and  vzdues  for  the  down¬ 
stream  droplets  are  profoundly  less  than  those  for  the  lead  droplet.  The 
transport  rates  for  the  third  droplet  au'e  lower  than  those  for  the  second 
droplet,  except  when  the  lead  droplet’s  interaction  with  the  second  droplet 
becomes  very  strong.  The  results  of  Case  4  show  that  the  second  and  third 
droplets  exhibit  almost  identical  behaviors.  The  interaction  between  the 
downstream  droplet  and  the  lead  droplet  becomes  very  significant  for  the 
small  initial  spacing  case. 

The  variations  of  trajectory  with  time  and  drag  coefficients  vs.  instan¬ 
taneous  Reynolds  number  for  the  cases  of  different  initial  spacings  are  illus¬ 
trated  in  Figure  8  (Cases  1,  2,  and  3)  and  Figure  9  (Cases  4,  5,  and  6), 
respectively.  Results  of  Case  1  with  large  initial  spacing  (D12=D23=K) 
indicate  that  the  drag  coefficient  is  smaller  for  the  second  droplet  than  for 
the  lead  droplet  and  still  smaller  for  the  third  droplet.  The  major  decrease 
in  drag  occxurs  in  the  first  two  droplets.  In  Case  2,  the  drag  coefficients  of 
the  lead  and  the  second  droplets  are  significantly  reduced  (comparing  curves 
7  and  10,  and  8  and  11,  respectively)  due  to  the  strong  interactions  when 
the  first  two  droplets  are  spaced  only  two  diameters  away.  Similar  trends 
occur  for  the  downstream  droplet  pairs  of  Case  3  (comparing  curves  8  and 
14,2uid  9  and  15,  respectively).  The  third  droplet  in  Case  2  has  a  higher 
drag  coefficient  than  that  of  the  second  droplet  since  the  second  droplet  is 
better  shielded  by  the  droplet  before  it.  The  D23  thus  increases  with  time. 
Also,  note  that  the  drag  coefficient  of  the  third  droplet,  which  is  spaced  far 
away  from  the  second  droplet,  seems  to  be  independent  of  D12  as  indicated 
in  curves  9  and  12.  However,  it  strongly  depends  upon  D23  as  illustrated  in 
curves  9,  12  and  15. 

The  comparison  between  results  of  Figures  8  and  9  leads  to  the  clear 
conclusion  that  as  initial  spacings  decrease,  the  transport  rates  decrease 
correspondingly.  The  drag  coefficients  for  the  downstream  droplets  are  2dl 
collapsed  together.  However,  at  the  final  calculations  when  three  droplets 
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interact  extensively  with  neighbors,  the  second  droplet  experiences  not  only 
downstream  interactions  from  the  lead  droplet  but  jJso  upstream  interaurtions 
from  the  third  droplet.  As  a  result,  the  second  droplet  possesses  the  lowest 
drag  coefficient.  Note  the  interaction  from  the  upstream  droplet  is  consider¬ 
ably  larger  than  that  from  the  downstream  droplet.  The  lead  droplet  in  Case 
2  has  the  lowest  drag  coefficient  (curve  10)  among  the  lead  droplets  of  three 
cases  since  the  interaction  from  the  second  droplet  is  the  strongest  (sp2u;ing 
D12  is  the  smallest).  The  behavior  of  the  spacing  variation  is  qualitatively 
similar  to  that  of  large  initial  spacing  case.  The  variation  of  D23,  for  the 
cases  of  approaching  droplets,  is  always  negligible  except  during  the  final 
calculation  period. 

In  order  to  investigate  the  effect  of  D23  on  the  first  and  the  second 
droplets  with  a  given  D12,  the  comparisons  of  drag  coefficients  for  the  first 
two  droplets  with  variable  D23s  are  presented  in  Figures  10  and  11.  The 
comparisons  of  results  with  the  same  D12  in  the  two-droplet-only  arrange¬ 
ment  are  also  presented.  For  the  cases  of  small  initial-Dl2,  the  main  droplet 
interactions  occur  between  the  first  two  droplets.  The  presence  of  the  third 
droplet  does  not  significantly  affect  the  lead  droplet;  it  reduces  the  Sherwood 
and  Nusselt  numbers  of  the  lead  droplet  by  less  than  5  %  in  magnitude.  The 
drag  coefficients  and  Nusselt  numbers  for  the  second  droplet  decrease  as  the 
D23  decreases,  though  the  percent  change  is  small.  For  the  cases  of  large  ini¬ 
tial  D12,  the  behavior  of  the  first  droplet  becomes  insensitive  to  the  presence 
of  the  third  droplet.  An  interesting  observation  is  that  the  presence  of  the 
third  droplet  at  a  sufficient  distance  from  the  second  droplet,  with  insignif¬ 
icant  hydrodynamic  interaction  between  the  second  and  the  third  droplet, 
may  increase  the  drag  coefficient  of  the  second  droplet  (compeuiing  curves  4 
and  5  in  Figures  11).  The  second  droplet  in  the  two-droplet-only  arrange¬ 
ment  receives  the  interaction  of  recirculating-thermal-entrainment  from  the 
free  stream;  as  a  result,  the  transfer  number  is  high.  Therefore,  the  surface 
blowing  is  enhanced  2md  the  friction  drag  is  reduced.  The  addition  of  the 
third  droplet  to  the  system  moves  this  thermal  effect  to  the  third  droplet; 
hence,  the  transport  rates  of  the  second  droplet  aie  recovered.  However,  for 
the  case  of  a  small  D23,  the  strong  interaction  from  the  third  droplet  reduces 
the  transport  rates  of  the  second  droplet  by  a  large  magnitude  (curve  6). 

A  nonlinear  regression  model  using  least  squares  has  been  employed  to 
find  the  correlations  between  the  transport  rates  of  an  interacting  droplet  auid 
the  corresponding  transport  rates  when  droplet  is  isolated.  The  generalized 
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form  of  numerical  correlations  for  drag  coefficient  and  Nusselt  number  2uid 
Sherwood  number  of  each  droplet  can  be  expressed  as 

For  the  lead  droplet: 

=  0.759Dl2°-*£>23°“^ 

^Di»o 

=  0.752i)12°““I>23°“'^ 

=  0.231D12-°=*“Z)23-°-2“ 

Sfl  filmic 

For  the  second  droplet: 

=  0.211i}12“-^Z)23°'^ 

^Diao 

-  o.296£>12°-^“Z>23°-*®® 

=  4.836I>12-°-‘®®i}23-®-^" 

Sh/ilrni^ 

For  the  third  droplet: 

=  0.434Z?12-°'“Z?23°“® 

Coiao 

-  o.343D12°  *“Z?23“-*“ 

^^Jilmiao 

_  io.476P12'°®®®I>23~°-^ 

Shfilmiao 

The  correlations  are  valid  for  droplet  spacing  ranged  from  2.7  to  12  and 
for  Reynolds  number  ranged  from  90  to  130.  Some  correction  factors,  with 
different  combinations  of  D  12s  and  D23s,  computed  from  above  correlations 
are  presented  in  Table  2.  The  correlations  basically  predict  the  right  trends 
of  interacting  effects  on  droplet  transport  rates.  The  influence  of  D 12  on  the 
third  droplet  is  insignificant  as  indicated  by  the  small  exponents  of  D12. 
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3  Conclusions 


The  detailed  behaviors  of  three  vaporizing,  interacting  droplets  for  the  cases 
of  different  initi<il  spacings  have  been  investigated  carefully.  The  results 
indicate  that  the  interacting  effects  are  strongly  dependent  upon  the  initial 
droplet  spacings.  The  general  qualitative  conclusions  drawn  from  the  two- 
droplet  study  [2]  can  be  applied  to  two  neighboring  droplets  in  the  three- 
droplet  cinadysis. 

Results  for  the  cases  of  sufficiently  large  spacing  (above  approximately  6 
droplet  diaimeters)  show  that  the  flow  field  of  each  droplet  is  qualitatively 
similar  to  that  of  an  isolated  droplet,  although  the  transport  rates  axe  reduced 
along  the  downstream  direction  due  to  the  cascade  effect  of  the  wake.  The 
major  drops  in  transport  rates  occur  in  the  first  two  droplets.  For  the  cases  of 
small  droplet  spacings  (less  than  approximately  3  droplet  diauneters),  the  flow 
field  of  downstream  droplets  can  be  significantly  altered  due  to  the  interzuztion 
effects  from  the  upstream  droplets.  Usually,  the  second  droplet  has  the  lowest 
drag  coefficient  since  it  receives  interactions  from  both  neighboring  droplets. 
However,  the  difiFerence  in  transport  rates  between  the  second  and  the  third 
droplets  is  not  significant  since  both  droplets  are  fully  protected  by  the  wake 
of  the  first  droplet.  The  effect  of  D23  on  the  behavior  of  the  first  droplet  is 
insignificant.  However,  depending  upon  the  values  of  D12  as  well  as  D23,  the 
effect  of  D23  on  the  behavior  of  the  second  droplet  may  become  significant. 
The  correlations  for  heat  transfer  and  droplet  dynamics  have  been  developed 
and  can  be  applied  in  a  droplet  spray  calculation. 
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Table  1:  Initicd  Droplet  Spacings  Considered  In  Each  Case  of  the  Three- 
Droplet  Study  _ 


Case 

D12 

D23 

1. 

12 

12 

2. 

4 

12 

3* 

12 

4 

4. 

6 

6 

5. 

4 

6 

6. 

6 

4 

Table  2:  Correction  Factors  for  Cd,  Nu  and  Sh  with  Different  Combinations 
of  D12  and  D23 


D12 

D23 

msnm 

Wfsrrm 

KS9 

1^3 

wraswm 

mmmi 

1.004 

0.781 

0.640 

0.874 

0.752 

0.687 

0.574 

1.001 

0.694 

0.540 

0.870 

0.657 

0.601 

RSil 

0.688 

0.872 

1.000 

0.647 

0.486 

0,868 

0.559 

1.077 

12.00 

0.935 

0.611 

0.659 

0.847 

0.663 

0.667 

0.738 

0.693 

4uM 

6.00 

0.931 

0.542 

0.563 

0.841 

0.579 

0.577 

0.927 

0.972 

1.137 

6.00 

4.00 

0.929 

0.510 

0.838 

0.535 

0.533 

0.918 

1.142 

1.465 

4.00 

12.00 

0.897 

0.529 

0.667 

0.832 

0.616 

0.659 

0.945 

0.749 

4.00 

6.00 

0.892 

0.470 

0.574 

0.825 

0.539 

0.566 

0.922 

1.189 

1.282 

4.00 

4.00 

0.890 

0.438 

0.522 

0.822 

0.498 

0.521 

0.911 

1.397 

1.696 

2.70 

12.00 

0.861 

0.460 

0.673 

0.819 

0.574 

0.653 

0.944 

1.099 

0.792 

2.70 

6.00 

0.857 

0.408 

0.582 

0.811 

0.502 

0.558 

0.918 

1.447 

1.404 

2.70 

4.00 

0.855 

0.381 

0.532 

0.808 

0.464 

0.512 

0.906 

1.700 

1.905 
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Figure  1:  Flow  field  configuration  and  grid  distribution  at  two  different  times. 
Figure  2:  Vorticity  and  isotherm  contour  plot  of  gas  and  liquid  phases  with 
different  initial  droplet  spacings. 

Figure  3:  Mass-fraction  contour  of  case  1  at  two  different  times. 

Figure  4:  Local  Nusselt  number  distribution  of  the  three  droplets  with  dif¬ 
ferent  initial  droplet  spacings. 

-  1  Lead  droplet,  Case  1 ,  Re(t)  =  85.6,  s.p.  at  143.2,  local  min.  at  1493 

- 2  Lead  droplet.  Case  4,  Re(t)  =  86.5,  s.p.  at  1303,  local  min.  at  141 A 

- 3  Second  droplet.  Case  1,  Re(t)  =  88.7,  s.p.  at  145.9,  local  min.  at  150.7 

-  4  Second  droplet.  Case  4,  Re(t)  =  913,  s.p.  at  156.6,  local  min.  at  0.0 

-  5  Third  droplet.  Case  1,Re(t)  =  89.5,  s.p.  at  150.1,  local  min.  at  158.8 

- 6  Third  droplet.  Case  4,  Re{t)  =  91.3,  s.p.  at  167.1,  local  min.  at  0.0 

Figure  5:  Surface  shear  stress  distribution  of  the  three  droplets  with  different 
initial  droplet  spacings. 

_  1  _Lead  droplet.  Case  1,  Re(t)  s  85.6,  s.p.  at  1433,  local  min.  at  121.7 

_  2  Lead  droplet.  Case  4,  Re(t)  =  86.5,  s.p.  at  1303,  local  min.  at  114.6 

_ 3  Second  droplet,  Case  1,  Reft)  s  88.7,  s.p.  at  145.9,  local  min.  at  1253 

_  4  Second  droplet,  Case  4,  Reft)  =  913,  s.p.  at  156.6,  local  min.  at  1493 

_  5  Third  droplet,  Case  1,  Reft)  =  89.5,  s.p.  at  150.1,  local  min.  at  129.4 

_ 6  Third  droplet.  Case  4,  Reft)  =  91.3,  s.p.  at  167.1,  local  min.  at  157.7 

Figure  6:  Time  variation  of  drag  coeflBcients  of  the  three  droplets  with  dif¬ 
ferent  initial  droplet  spacings. 

_  1.  Lead  Droplet,  flnftlal  012=12,023=12) 

_ 2-  Second  Droplet,  finitlal  012=12,  023=12) 

_ 3.  Third  Droplet,  flnltlal  012=12,023=12) 

_  4.  Lead  Droplet,  flnltlal  012=6,023=6) 

_  5.  Second  Droplet,  flnltlal  012=6,  023=6) 

_ 6.  Third  Droplet,  flnltlal  012=6,  023=6) 

_ 7.  Isolated  Droplet 


Figure  7;  Time  vaxiation  of  Nusselt  numbers  of  the  three  droplets  with  dif¬ 
ferent  initial  droplet  spacings. 

-  1.  Lead  Droplet,  (Initial  Di2ai2,  D23=12) 

- 2.  Second  Droplet,  (Initial  D12=12,  D23=12) 

- 3.  Third  Droplet,  (Initial  Di2si2,  D23sl2) 

-  4.  Lead  Droplet,  (Initial  D12=6, 023=6) 

_  5.  Second  Droplet,  (Initial  D12=6,  D23=6) 

- 6.  Third  Droplet,  (Initial  D12=6,  D23=6) 

_ 7.  Isolated  Droplet 

Figure  8:  Time  variation  of  droplet  drag  coefficients  and  droplet  spacings  for 
the  cases  of  large  initial  droplet  spacings  (Cases  1,  2  and  3) 

1.  D12,  Case  1;  2.  D23,  Case  1;  3.  Dl2,  Case  2; 

4.  D23,  Case  2;  5.  D12,  Case  3;  6.  D23,  Case  3; 

7.  Cn  for  the  Lead  Droplet,  Case  1;  8.  Cj,  for  the  Second  Droplet,  Case  1: 

9.  Cd  for  the  Third  Droplet,  Case  1;  10.  Co  for  the  Lead  Droplet,  Case  2- 

Jr  n  Droplet,  Case  2;  12.  Cj,  for  the  Third  Droplet,  Case  2; 

13.  Co  for  the  Lead  Droplet,  Case  3;  14.  Co  for  the  Second  Droplet,  Case  3; 
15.  Co  for  the  Third  Droplet,  Case  3;  16.  Co  for  an  Isolated  Droplet; 


Figure  9;  Time  v2Lriation  of  droplet  drag  coefficients  and  droplet  spacings  for 
the  cases  of  large  initial  droplet  spacings  (Cases  4,  5  and  6) 

I.  D12,  Case  4;  2,  D23,  Case  4;  3.  D12,  Case  5; 

4.  D23,  Case  5;  5.  Dl2,  Case  6;  6.  D23,  Case  6; 

7.  Co  for  the  Lead  Droplet,  Case  4;  8.  Co  for  the  Second  Droplet,  Case  4; 

9.  Co  for  the  Third  Droplet,  Case  4;  10.  Co  for  the  Lead  Droplet,  Case  5; 

II.  Co  for  the  Second  Droplet,  Case  5;  12.  Co  for  the  Third  Droplet,  Case  5; 
13.  Co  for  the  Lead  Droplet,  Case  6;  14.  Co  for  the  Second  Droplet,  Case  6; 
15.  Co  for  the  Third  Droplet,  Case  6;  16.  Co  for  an  Isolated  Droplet; 


Figure  10:  Time  vaxiatioa  of  drag  coefBcients  of  the  first  two  droplets  for 
D12  =  4  and  variable  D23s 

_  1.  The  First  Droplet  of  Two-Oroplet  Calculation,  D12=4 

_ 2-  The  Rrst  Droplet  of  Three-Droplet  Calculation,  012=4,  D23=12 

_ 3.  The  Rrst  Droplet  of  Three-Droplet  Calculation,  D12=4,  D23=6 

_  4.  The  Second  Droplet  of  Two-Droplet  Calculation,  Dl2s4 

_  5.  The  Second  Droplet  of  Three-Droplet  Calculation,  D12=4,  D23=12 

_ 6.  The  Second  Droplet  of  Three-Droplet  Calculation,  D12=4,  D23=6 

Figure  11:  Time  variation  of  drag  coefficients  of  the  first  two  droplets  for 
D12  =  12  and  vziriable  D23s 

_  1.  The  Rrst  Droplet  of  Two-Droplet  Calculation,  D12=12 

_ 2.  The  Rrst  Droplet  of  Three-Droplet  Calculation,  D12=12,  D23=12 

_ 3.  The  First  Droplet  of  Three-Droplet  Calculation,  D12=12,  D23=4 

_  4.  The  Second  Droplet  of  Two-Droplet  Calculation,  D12=12 

_  5.  The  Second  Droplet  of  Three-Droplet  Calculation,  D12=12,  D23=12 

_ 6.  The  Second  Droplet  of  Three-Droplet  Calculation,  D12=12,  D23s4 


FLOW  CONFIGURATION  -  SCHEMATIC 


Time  = 

Rx  =  1.000,  iZj  =  1.000, /23  =  1.000, 

D12=6.00,  D23=6.00,i?ei=i2e2=i?e3=100. 


D12=2.68,  D23=6.01,  i?e, =84.33,  /?e2=90.26, 


8.8  15.7  22.7  29.6  36.6 

Time  =  3.00,  Re^  =  79.93,  Rci  =  84.50,  =  85.90 

Rx  =  0.998,  iZj  =  0.999,  i?3  =0.999,  D12=8.73,  D23=11.14 

Initial  R^  =  1.0,  R^  =  1.0,  .R3  =  1-0,  D12=12,  D23=12 


8.8  14.0  19.2  24.4  29.6 

Time  =  5.00,  Re^  =  70.55, i2e2  =  77.52,  Re^  =  79.11 
Ri  =  0.991,  i?2  =  0.995,  iZa  =  0.996,  Dl2=2.96,  D23=9.50 
Initial  Ri  =  1.0,  i?2  =  1.0,i23  =  1.0,  D12=12,  D23=12 


8.9  13.2  17.4  21.7  26.0 

Time  =  2.00,  Ra  =  86.53,  Rej  =  91.29,  Res  =  91.28 
Ri  =  0.998,  .Rj  =  0.999,^3  =0.999,  D12=3.74,  D23=5.93 
Initial  Ri  =  1.0,  R2  =  1.0,R3  =  1.0,  D12=6,  D23=6 


</)  -0.88 


0.  18.  36.  54.  72.  90.  108.  126.  144.  162.  180. 
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Abstract 

The  detailed  analysis  of  a  cold  LOX  droplet  suddenly 
introduced  into  the  hot  methane  fuel-vapor  stream  has 
been  performed  in  this  study.  The  effects  of  variable 
thermophysical  properties,  real  gas  behavior,  transient 
heating  eind  internal  circulation  of  liquid,  deceleration 
of  the  flow  due  to  the  drag  of  the  droplet,  boundary- 
layer  blowing,  and  moving  interface  are  included.  A 
primitive-variable  formulation  with  an  implicit  finite- 
difference  scheme  has  been  developed  to  solve  the  com¬ 
plete  set  of  Navier-Stokes,  energy,  and  species  equa¬ 
tions.  The  results  are  presented  for  the  low  pressure 
case.  The  interesting  phenomena  due  to  the  large  sur¬ 
face  blowing  and  surface  boiling  is  examined.  The  mod¬ 
ifications  required  for  the  current  model  to  treat  high 
pressure  case  are  discussed. 
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Greek 

e 

<x 

Kg 

K\ 

Ps 

Pi 


t'  ja!^,  radial  coordinate 
R’‘T^I{U'^  q  M'^),  gas  constant 
universal  gas  constant 

«oU«.o/^,oo//^i,=o.  gas-phase  Reynolds  number 
“qUto.oP/  o/P/  o>  liquid-phase  Reynolds  number 
T'/T^,  temperature 
time 

Y'lU'aa  Q,  velocity 

mass  fraction 

z' ja'^,  axial  coordinate 


angular  coordinate 

Qa'o),  surface  tension 
Kg/Kg  oo,  conductivity  of  gas  phase 
0,  conductivity  of  liquid  phase 
P'glP's.oo’  viscosity  of  gas  phase 
p't/p'i  o>  viscosity  of  liquid  phase 
P'j/Pi,oo>  density  of  gas  phase 


Cp, 

^p'gl^p'g  oo>  specific  heat  of  gas  phase 
Cp'JCp'i  g,  specific  heat  of  liquid  phase 

TEg 

^'l^'g.oo/Wo  Pj,ooQ>j,oo)’ 
gas-thermal-diffusion  time 

Vg 

F 

'^'gl'^'g,oo’  niass  diffusivity  of  gas  phase 
^7(Ui,,o^ao^p',,c»).  drag  force 

'I’E! 

P;,o^p/,o)> 

liquid-thermal-diffusion  time 

f 

h 

,  fugacity 

^'li<^p'g.oo'^g.oo)’  enthalpy 

I'Hg 

<X.=o/KVa,oo). 
gas-hydrodynamic-diffusion  time 

L 

LC  g 

L' f{Cp'i  qT^),  latent  heat  of  vaporization 
P',,oo^i,ooCpj,oo/«J,oo- 

thi 

Pl,o)i 

liquid- hydrodynamic-diffusion  time 

gas-phase  Lewis  number 
M  M  /M^o ;  M'  =  ’ 

equivalent  molecular  weight 
Pg  (p'g  -  P'g.co)/iPg,coU!„/),  gas-phasc  pressure 

PI  (P;  -  P't,o)/(p'i.o^L/)’  liquid-phase  pressure 

Peg  RegPrg,  gas-phase  Peclet  number 

Pei  ReiPri,  liquid-phase  Peclet  number 

Pi'g  P'g,ooCp'g  ,^lK'g  ^,  gas-phase  Prandtl  number 

Prt  P/,o^p;  o/'*/,0'  liquid-phase  Prandtl  number 
Q'  heat  flux 
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Ts  ^/(oq  ),  gas-species-diffusion  time 

<i>  fugacity  coefficient 


Subscripts 

ave  volumetric  average 
crii  critical  point 

/  fuel 

g  gas  phase 

/  liquid  phase 

n  normal  direction 

o  oxygen 

si  saturated  liquid 
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V  vapor 

0  initial  conditions 

6  tangential  direction 

00  free  stream  conditions 

Superscript 

0  ideal  gas  state 

'  dimensional  quantity 

1  Introduction 

The  study  of  the  liquid  behavior  of  liquid-propellant 
rocket  engines,  where  both  liquid  fuel  and  oxidizer  are 
injected  from  one  end  of  the  chamber,  has  recently  be¬ 
come  the  subject  of  various  experimental  and  theoret¬ 
ical  investigations.  One  of  the  characteristics  whereby 
liquid-fueled  rocket  engines  differ  from  air-breathing 
engines  is  that  both  fuel  droplets  and  liquid  oxygen 
(LOX)  droplets  are  present.  The  use  of  LOX/hydrogen 
fuel  propellant  combination  or  LOX/methane  propel¬ 
lant  combination  in  an  advanced  launch  vehicle  booster 
engine  appears  extremely  attractive  due  to  high  propel¬ 
lant  bulk  density  and  the  relatively  high  performance 
characteristics  of  these  propellant  combinations.  Hy¬ 
drogen  and  methane  are  cryogenic  fuels  and  would 
be  injected  in  a  fluid  state.  Usually  in  the  practical 
combustor,  liquid  hydrogen  or  liquid  methane  vaporize 
faster  than  the  oxygen  droplets.  As  a  result,  oxygen 
droplet  vaporization  becomes  a  rate-controlling  factor 
in  determining  the  mixing  performance  of  the  hetero¬ 
geneous  fuel/oxidizer  mixture. 

Due  to  some  technical  difficulties  in  handling  cryo¬ 
genic  fuels  €Uid  LOX  droplets  at  extremely  low  tem¬ 
perature  in  the  regular  laboratory,  only  very  lim¬ 
ited  experimental  works  [1,  2]  utilized  'LOX/H2  and 
L0X/C/f4  propellant  to  study  combustion  perfor¬ 
mance  have  been  constructed.  Although  substantial 
theoretical/computational  results  have  been  reported, 
the  theoretical  studies  of  the  vaporization  of  LOX 
droplets  have  not  reached  a  satisfactory  status.  Many 
assumptions  employed  to  simplify  analysis  have  ren¬ 
dered  the  results  to  become  questionable  when  detailed 
transport  processes  of  LOX  vaporization  occur  at  a 
high  temperature,  high  pressure  and  convective  envi¬ 
ronment. 

Generally,  the  strategy  of  droplet  research  over  the 
past  decade  has  been  to  develop  simplified  transient 
vaporization  models  which  are  physically  accurate  but 
sufficiently  computationally  efficient  so  that  hundreds 
or  thousands  of  droplets  in  a  spray  combustor  can  be 
tracked  [3,  4,  5].  During  the  past  decade,  the  research 
of  LOX  vaporization  has  more  or  less  followed  the  same 
trend.  The  existing  LOX  models  are  mainly  focused  on 
the  one-dimensional  calculations  where  quasi-steady  as¬ 
sumption  and  film  theory  are  employed  for  gas-phase 


analysis  while  simplified  conduction  models(  such  as  in¬ 
finite  conductivity,  conduction  limit,  and  effective  con¬ 
ductivity  models)  are  frequently  used  for  liquid-phase 
analysis.  No  convection  and  variable  properties  effects 
have  ever  been  considered  in  the  existing  LO.X  mod¬ 
els.  It  is  well  known  that  the  quasi-steady  assumption 
becomes  invalid  with  the  high  pressure  environment. 
The  internal  motion  of  the  liquid  phase,  induced  by  the 
gas-phase  convection,  may  significantly  change  the  time 
scales  of  transport  processes  of  the  droplet.  Variable 
property  effects  due  to  the  variations  of  temperature 
and  species  composition  can  also  affect  the  computed 
quantities.  These  effects  are  no  longer  negligible  in  the 
development  of  a  comprehensive  numerical  model  for 
the  LOX  droplet. 

In  parallel  with  the  development  of  simplified  LOX 
vaporization  modeling,  there  is  a  need  to  pursue  “e.x- 
act”  solutions  (via  finite-difference  calculations)  of  the 
flow  and  thermal  fields  surrounding  and  within  vaporiz¬ 
ing  LOX  droplets.  These  Navier-Stokes  solutions  serve 
two  major  purposes:  1)  correlations  for  quantities  such 
as  droplet  drag  coefficients  are  obtained  that  can  be 
used  in  spray  combustion  analyses  with  the  simplified 
models  (which  do  not  independently  predict  drag  co¬ 
efficients)  and  2)  the  exact  solutions  can  be  used  as  a 
standard  for  comparison  of  the  simplified  models.  Also 
the  numerical  data  can  be  provided  as  a  benchmark  for 
the  experimental  research. 

The  present  research  involves  the  <letailed  investi¬ 
gation  of  an  isolated  liquid  oxygen  droplet  vaporizing 
in  the  high  temperature  fuel-vapor  environment.  The 
forced  convection  of  the  gas  phase,  the  transient  decel¬ 
eration  of  the  flow  due  to  the  drag  force,  the  surface  re¬ 
gression,  the  internal  circulation  and  transient  heating 
of  the  liquid  phase,  and  variable  properties  are  consid¬ 
ered.  The  LOX  vaporization  at  near  critical  conditions 
will  be  addressed  since  most  rocket  engine  combustors 
are  operated  in  the  very  high  pressure  domain  where 
the  state  equation  and  thermodynamic  functions  are 
considerably  different. 

Our  current  axisymmetric  models  have  been  devel¬ 
oped  for  hydrocarbon  fuel  droplets  vaporizing  under 
subcritical  conditions  [6].  Therefore,  the  well  devel¬ 
oped  isolated-droplet-code  will  be  utilized  as  the  base 
code  to  facilitate  the  present  research.  The  problem 
features  very  complex  property  calculation  along  with 
very  high  transfer  number  (  «  0(100)  )and  high  density 
ratio.  The  early  rapid  surface  heating  is  worthwhile  to 
develop  a  new  scheme  to  surmount  the  numerical  diffi¬ 
culties.  In  the  present  research,  some  effort  is  required 
to  adjust  these  transient  droplet  heating  and  vaporiza- 
tion  models  so  that  liquid-propellant  rocket  combustion 
can  be  studied. 

The  variation  of  liquid-phase  density  with  tempera¬ 
ture  is  significant.  The  LOX  may  experience  the  ther¬ 
mal  expansion  during  the  transient  heating  and  the 
surface  moving  velocity  cannot  be  neglected.  In  or- 
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der  to  carefully  account  for  unsteadiness  due  to  density 
variation,  the  primitive  variables  (K,  and  pi)  for¬ 
mulation  for  the  liquid  phase  is  employed.  Note  that 
in  previous  codes  [6,  7],  the  liquid-phase  momentum 
equations  are  formulated  by  a  streamfunction-vorticity 
approach.  The  modifications  of  the  code  to  adapt  the 
primitive  variables  for  the  liquid  phase  as  well  as  to 
employ  more  accurate  interface  conditions  have  been 
developed.  The  details  will  be  discussed  in  the  follow¬ 
ing  sections. 

2  Formulation 

The  schematic  flow  configuration  is  presented  in  Fig¬ 
ure  1  where  the  flow  passing  over  a  vaporizing  LOX 
droplet  moving  with  respect  to  gas  stream  is  shown. 
The  flow  is  laminar  and  axisymmetric  with  initially  uni¬ 
form  ambient  conditions  specified  by 
and  Vo  00  =  0.  The  frame  of  reference  is  fixed  to  the 
center  of  the  LOX  droplet.  The  problem  can  be  viewed 
as  an  impulsively-started  flow  over  a  fixed  LOX  droplet. 
Unless  otherwise  stated,  the  basic  assumptions  remain 
the  same  as  mentioned  in  Chiang  et  al  [6]. 

The  conservation  equations  in  both  phases  are  writ¬ 
ten  in  non-dimensional  forms  with  respect  to  a  cylin¬ 
drical  coordinate  system.  The  initial  radius,  upstream 
velocity  and  thermophysical  properties  have  been  used 
to  non-dimensionalize  the  variables.  The  diffusion 
time  has  been  selected  as  the  time  scale  in  this  study. 
According  to  the  non-dimensionalization  given  in  the 
nomenclature,  the  governing  equations  are  presented 
below. 

2.1  Governing  Equations 


Gas  Phase 


Continuity  Equation 


+  y^{rRe,p,Vr)  +  ^JrRe,p,V,)  =  0  (I) 
Momentum  Equation  in  r-direction 


2d.  ..dK  Vr  dV,..d.  rdVr^dV, 


+  RcgPg  - 


3  ^  r  dr 


(2) 


Momentum  Equation  in  z-direction 

I - {rp,V,)  +  ^(rRegPgVrV,) 

OTHg  or 


+^JrRegipgV.V.+Pg)}  = 


,  r^dV,  Vr  dVr 

-  --  -  +  Fd  (3) 


where  Fq  is  the  D’Alembert  force  (reversed  inertial 
force  due  to  the  drag  force  on  the  droplet)  which  is 
uniformly  applied  to  the  whole  gas-phase  flowfield. 
Energy  Equation 


+(/»<,  -  /»/) 


^^irp9yo)  +  ^irPe,p,Y.Vr) 


+^^irPegPgY,V,) 


d,  dTg.d.  dTg. 


~{rLegiC,^-C,JTgPgVg^) 


-  yjUg{C,j  -  CMp.Vg^)  (4) 

where  h  =  Cp^dT 
Species  Equation 


drs 


dz^  Leo 


d  .  .r,dYi.  d  .  _  dYi. 


ar'" dr 
Equation  of  State 


(5) 


(6) 

We  assume  that  the  mixture  behaves  as  an  ideal 
gas.  This  is  a  good  assumption  for  the  low  pressure 
case.  The  assumption  may  become  impractical  for 
pressure  near  the  critical  point.  The  other  alterna¬ 
tive  to  relate  Pg,P},Tg  and  Vi  is  to  employ  the  two- 
parameter  Redlich-Kwong  equation  [8]  with  the  mixing 
rule  of  Chueh  and  Prausnitz  [9].  Our  test  runs  indi¬ 
cated  that  the  computational  time  will  increase  tremen¬ 
dously  since  it  is  necessary  to  solve  the  real  root  of  a 
third-order  polynominal  equation  which  is  a  very  time- 
consuming  task  within  an  iterative  Poisson  equation 
solver.  For  the  time  being,  in  order  to  prevent  e.xhaust- 
ing  our  precious  supercomputer  allocation,  we  employ 
this  assumption  for  the  pressure  below  40  atms. 
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Liquid  Phase 


Continuity  Equation 

+  §^{rRe,p,Vr)  +  ^SrRtipiV,)  =  0  (7) 
Momentum  Equation  in  r-direction 
^—^{rpiVr)+-^{rRei{p,VrVr+pi)}  +  ^{rRe,piVrVi) 

2d.  ..dVr  Vr  dK..d.  rdVr^dV,.. 

+  2/i,  v;  dVr  dV, 

+  Re,p,-—(2y--^--^)  (8) 

Momentum  Equation  in  z-direction 
^(rp,V,)+^(rRe,p,VrK)  +  l{rRe,(p,V.V,+p,)} 


d  ,  ,dVr  dV;„  2  d  ,  ,^dV,  Vr  dVr,, 

=  ^  +  3  -  T  - 


(9) 


Cp,  { +  ^irPe,p,VrT,)  +  ±{rPe,p,VTi)} 


Energy  Equation 


d  ,  dT,.  d  .  dTu 


(10) 


Pressure- Density-Temperature  Relation 
The  Hankinson-Brobst-Thomson  technique  is  used  to 
predict  the  compressed  liquid  density. 

‘  ■  iiW]  <“> 


^  P'1., 


1  —  cln 


where  p/,a,,  the  saturated  liquid  density  at  the  vapor 
pressure,  j/^p.  j3  is  obtained  from 

0/p'crit  =  -1  +  a(l  -  +  6(1  -  r.)2/3  +  d(l  -  Tr) 

+  e{l-Tr)‘'^^  (12) 

where  Tr  =  Values  of  constants  are  available  in 

^  erf« 

Reid  et  al  [10].  We  assume  that  there  is  no  absorption 
of  fuel- vapor  into  the  liquid  phase;  as  a  result,  there  is 
no  need  to  solve  species  equations  for  the  liquid  phase. 

2.2  Initial  and  Boundary  Conditions 


2.2.1  Initial  Conditions 

Gas  Phase:  Vr  -Pg  -Yg  =  0,Vi  =Tg  =  Pg  =■  I 
Liquid  Phase;  Vr  =  pi  =  V^  =  0,pi  =  1,T;  =  T(,o 
Gas/Liquid  Interface  :  Vr  =  Vz  =  Yo  =  Pg  =  pi  = 
Q,Tg  -Ti  —  Tio,Pj  =  ^/Ti^o.pt  =  1 


2.2.2  Gas/Liquid  Interface  Boundary  Condi¬ 
tions 

Tangential  Stress  Condition: 

Ke.g  —  '’ili.l  =  ^ 

where  cr'  =  <t[,  -f-  ^T/;  Usually  ^  is  assumed  to  be  a 
constant.  The  nondimensional  form  becomes 


JtL 

Rei 


(dVl,g  Vi,e  ldV,,n\ 

V  dn  a  a  dO  ) , 


Pj,oo  Pg 


(dVg,f  Vg,g  I  dVg,„\  ^  I  d<r  dT, 
V  dn  a  a  de  J,  adT,  dO 


Pl.O  R^3 
Normal  Stress  Condition  : 

P’l  +  <n.l-p'g-  Kn.g  =  2<T7a' 


(13) 


The  corresponding  nondimensional-form  of  the  above 
equation  is 


p'g.oo  _  Ci  pi  dV;, 


^  Rt,  dn 


—  ^  2^'°°  dVy^n 


=  2--b 
a 


pj  Q  dfi 

^,oo  p{,o 
^oo.O^Pt.O 


Continuity  of  Tangential  Velocity: 

Yg.$.z  —  Vl^ff^z 


(14) 


(15) 


Conservation  of  Mass  Flux: 

By  assuming  no  accumulation  of  mass  at  surface,  the 
integr2il  form  of  the  continuity  equation  can  be  written 


as 


DU  f  D  u  ^  ^ 

RCgVg  „  —  —  (  RCgV,  „  — - 

P',.,  V  <^'^Hg 

where  the  droplet  regression  rate  is  give 
d  a 


(16) 


dTHi 


da  1 
dTHo  Re, 


sinOde  -t- 


(1- 


3  dVHg 

Note  that  this  equation  is  employed  to  calculate  the 
instantaneous  radius. 

Continuity  of  Temperature: 


Tg.i  =  T., 

Conservation  of  Energy: 


(18) 


\k\)  dn^>  dn'‘ 
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+  ‘  )].f^ 

Kl  \^,0 

Conservation  of  Species: 


dn),~  ^^9Scg[iVg,„  -  — — )(n  -  1)], 


Yo  —  T,  Relation  for  Phase  Equilibrium: 

The  vapor  pressure  is  usually  estimated  as  funi.Uon  of 
surface  temperature.  The  correlation  taken  from  Reid 
et  al.  [10]  is  employed  in  the  current  study. 

HPv.o/p'crit)  =  (1  -  xy'-lAx  +  Bx^-^  +<71^  +  £»x®] 

(21) 

where  x  =  1  - 

This  kind  of  correlation  have  been  particularly  suc¬ 
cessful  in  computational  programing. 

The  mass  fraction  at  the  interface  can  be  determined 
by 

At  low  pressure  case,  the  ideal-solution  assumption 
is  valid  and  the  Lewis  fugacity  rule  leads  to  the  Raoult’ 
law  which  simply  says  that 


However,  the  ideal  behavior  of  gas  mixture  no  longer 
occurs  under  the  high  pressure  situation.  The  va.- 
por/liquid  equilibrium  requires  that  fugacity  of  each 
species  component  must  be  equal  in  both  phases  at 
the  surface;  that  is  =  /,(,  where  subscript  i  stands 
for  species.  Hence,  it  is  necessary  to  compute  fugacity 
of  oxygen  in  both  phases.  The  mole  fraction  of  oxy¬ 
gen  vapor  can  be  determined  once  fugacity  is  available. 
Since  the  fugacity  has  a  dependency  on  mole  fraction  ( 
through  the  mixing  rule),  the  procedure  for  the  calcu¬ 
lation  of  Xo  requires  the  iteration  of  Xo  until  a  conver¬ 
gence  criteria  is  satisfied.  The  equations  for  iteration 
are  given  in  the  following  two  equations. 


Pv,o  ^st,o 

p's  <t>v.o 


v,dp 


where  Vi  is  the  molar  density  of  the  liquid  phase. 

Vapor-phase  fugacity  coefficient  of  oxygen,  4>v,o,  is 
given  by 

„5) 

In  order  to  integrate  this  equation,  the  Redlich- 
Kwong  equation  of  state  combined  with  the  mixing 
rules  of  Chueh  and  Prausnitz  must  be  employed.  Xo 
is  implicitly  involved  in  the  above  equation. 


The  enthalpy  of  vaporization  for  oxygen  can  be  cal¬ 
culated  from  the  thermodynamic  relationship 

_  d\n4>v^o  pg, 

Ro'T-  dT  ’ 

B  a  —  hi  a  (27) 

The  enthcilpy  of  vaporization  at  high  pressure  case 
is  expected  to  be  different  from  the  conventional  latent 
heat  which  is  a  function  of  temperature  only. 

The  detailed  computations  of  droplet  thermophysical 
properties  for  the  high  pressure  case  Ccin  be  found  in 
Reference  [12,  13,  14,  15,  16,  17] 

2.2.3  Outflow  (r  =  Too ,  t/2  <  0  <  w)  Boundary 
Conditions: 

^(v'.)  =  .g(v,)  =  :gm  =  ^(p,)=  °(y.)=o 


2.2.4  Inflow  (r  =  roo,0  <  d  <  'x/2)  Boundary 
Conditions: 

p,  =  n  =  V;  =  0,T  =  /),  =  14  =  1  (29) 

2.2.5  Axis  of  Symmetry  (0  <  r  <  Too,  ^  =  0,  t) 
Boundary  Conditions 


Gas  Phase: 


dVn  _dPg  _dT  _  dpg  _  dYc 


"  ~  89  ~  d9  ~  89  ~  89 
Liquid  Phase: 


^=0  (30) 


~  89  ~  89  ~  89  ~  89  ~  ^  ^ 

2.2.6  Boundary  Conditions  for  Governing 
Equations  of  Liquid  Phase 

At  Droplet  Surface 

Solutions  obtained  from  the  boundary-conditions- 
solver  are  regarded  as  the  Dirichlet  conditions  for  the 
liquid-phase  solver. 

At  Droplet  Center 

8Ti  _  8V.  _  8Vr  _dpi 
8n  dn~8n~8n~  ^  ’ 

3  Solution  Procedure 

A  brief  description  of  the  computational  methodologies 
is  given  here.  The  finite-difference  equations,  and  de¬ 
tailed  numerical  procedures  are  given  in  Chiang  [11]. 

Since  we  have  considered  internal  circulation  and 
transient  heating  of  the  liquid  phase,  forced  convection 
of  the  gas  phase,  transient  deceleration  of  the  flow  and 
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variable  properties,  it  is  then  necessary  to  solve  simul¬ 
taneously  the  complete  set  of  unsteady  Navier-Stokes 
equations,  energy  and  species  equations,  combined  with 
the  appropriate  boundary  conditions.  The  nonlinear 
and  highly  coupled  equations  make  the  analytical  solu¬ 
tion  almost  impossible.  We  have  to  resort  to  an  implicit 
finite-difference  numerical  algorithm. 

The  governing  equations  are  represented  in  general¬ 
ized  coordinates  which  conform  to  changing  boundaries 
due  to  decreasing  droplet  radii  associated  with  liquid 
component  evaporation.  Pressure  correction  equations 
are  employed  to  satisfy  indirectly  the  continuity  equa^ 
tions  in  both  phases.  The  pressure  correction  equa¬ 
tions,  which  are  Poisson  type  of  equations,  are  solved 
by  the  successive-over-relaocation  (SOR)  method.  Note 
that  the  equation  of  state  is  incorporated  with  the 
Poisson  equation  to  relate  pressure  and  density  of  gas 
phase.  The  density  is  immediately  updated  when  the 
new  estimate  of  pressure  becomes  available  during  the 
iterations.  In  the  pressure  correction  equation  of  liq¬ 
uid  phase,  we  temporarily  decouple  the  relationship 
between  density  and  pressure  to  avoid  the  undesired 
acoustic  pressure  wave.  Actually  the  liquid  density  is 
a  very  weak  function  of  pressure.  Hence,  the  varia^ 
tion  of  pressure  in  the  liquid  phase  is  assumed  to  be 
caused  by  the  fluid  motion  only.  The  density  is  updated 
until  the  iterations  of  pressure  correction  equation 
are  done.  The  momentum,  energy  and  species  equar 
tions  are  solved  by  an  alternating-direction-predictor- 
corrector  (ADPC)  method.  The  non-linear  gas/liquid 
interface  boundary  equations  are  treated  by  a  quasi- 
linearization  technique  and  solved  directly  by  the  in¬ 
version  of  tridiagonal  block  matrices.  The  governing 
equations  of  motion  as  well  as  the  interface  boundary 
conditions  are  solved  sequentially  in  an  interactive  se¬ 
quence  until  convergence  is  achieved  for  each  time-step 
of  the  calculation. 

Since  the  reference  frame  is  fixed  to  the  liquid 
droplet,  the  evaluation  of  the  drag  force  and  its  as¬ 
sociated  velocity  correction  for  the  gas  phase  are  incor¬ 
porated  in  the  iterative  process.  In  order  to  maintain 
a  dense  grid  distribution  at  the  droplet  interface,  the 
grid  locations  have  to  be  adjusted  at  each  time  step 
to  accommodate  droplet  surface  regression.  Obviously, 
the  metrics  of  the  transformation  have  to  be  updated 
whenever  the  grid  system  is  moved. 

In  the  overall  procedure,  the  sequential  solutions  of 
governing  equations  and  boundary  conditions  with  grid 
and  relative  velocity  adjustment  are  iterated  until  con¬ 
vergence  is  achieved.  After  convergence  is  reached,  the 
drag  coefficients,  average  Nusselt  and  Sherwood  num¬ 
bers  are  evaluated  at  prescribed  time  intervals. 

The  actual  property  computation  is  extremely  time- 
consuming  since  most  of  the  thermo  properties  are 
strongly  coupled  together  through  the  set  of  thermo¬ 
dynamics  equations.  An  iterative  approach  has  to  be 
employed.  We  have  compiled  the  available  property  in¬ 


formation  from  [10,  18,  19,  20]  and  constructed  the 
property  correlations  based  on  the  ambient  pressure 
and  an  extensive  range  of  temperature  before  start¬ 
ing  droplet  computation.  W’e  learned  that  some  of  the 
involved  material  properties  change  abruptly  near  the 
interface  where  the  temperature  gradient  is  large. 

The  major  numerical  difficulties  occur  at  a  very  early 
time.  The  large  gradients  at  interface  seem  to  create 
inappropriate  initial  profiles  used  to  start  iterations. 
Since  the  initial  conditions  are  physically  meaningless 
before  a  residence  time,  it  is  safe  to  employ  any  suitable 
profile  for  variables  in  order  to  form  a  converged  solu¬ 
tion.  We  have  designed  a  “ramp  temperature  profile,” 
where  the  ambient  temperature  increases  slowly  with 
time,  to  smooth  the  transient  gradient  of  temperature 
variation.  Also,  we  have  developed  a  control-volume 
scheme  to  solve  conservation  equations  at  the  interface 
sequentially  to  decouple  temporarily  variables  at  the 
beginning  of  computation.  After  a  near  converged  so¬ 
lution  is  obtained,  we  switch  to  the  tri-diagonal  block 
solver.  The  danger  of  solving  interface  conditions  se- 
quentiidly  is  that  the  temperature  may  shoot  above 
the  boiling  temperature  dramatically  when  the  LOX 
droplet  temperature  approaches  the  boiling  point. 

4  Results  and  Discussion 

The  numerical  simulations  of  an  isolated  LOX  droplet 
suddenly  injected  into  a  methane-fuel-vapor  flow  are 
studied  in  the  present  research.  The  ambient  temper¬ 
ature  and  initial  droplet  temperature,  are  selected  to 
be  1,000  k  and  100  k,  separately.  The  initial  droplet 
Reynolds  number  is  100.  The  critical  pressure  and  tem¬ 
perature  for  the  LOX  droplet  are  154.58  k  and  50.43 
aims,  respectively.  In  order  to  handle  the  surface  blow¬ 
ing  at  the  interface,  we  have  to  employ  a  very  small 
time  step  which  render  the  whole  computation  become 
very  time-consuming.  For  this  reason,  only  very  lim¬ 
ited  production  runs  with  pressure  ranging  from  10  to 
50  atms  are  made  in  order  to  study  the  pressure  ef¬ 
fects.  Unfortunately,  our  code  fails  when  the  ambient 
pressure  exceeds  40  atms.  The  detailed  reasons  will  be 
discussed  later.  We  only  present  results  from  the  low 
pressure  calculations. 

In  general,  the  flowfield  structure  is  similar  to  that 
of  a  fuel  droplet  vaporizing  in  an  air/fuel  mixture.  The 
detail  discussion  of  the  global  flowfield  is  given  in  Chi- 
ang  et  al  [6].  It  is  worthy  to  note  that  LOX  vapor¬ 
ization  features  highly  surface  blowing  and  highly  tran¬ 
sient  heating  which  leads  to  surface  boiling  at  very  early 
lifetime.  The  gas  flowfield  near  the  droplet  interface  is 
expected  to  be  significantly  influenced  by  the  vapor¬ 
ization.  The  transport  rates  are  also  expected  to  be 
significantly  reduced  due  to  the  surface  blowing. 

Several  snap  shots  of  the  flowfield  are  presented  in 
Figure  2  to  Figure  6.  Figure  2  shows  the  instan- 
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taneous  velocity  vectors  at  5  hydrodynamic  diffusion 
times  when  20%  of  mass  has  been  vaporized.  The  out¬ 
ward  velocity  vectors  are  observed  at  the  4  grid  points 
next  to  the  surface  of  the  front  stagnant  region.  Actu¬ 
ally,  the  surface  blowing  has  pushed  the  velocity  stag¬ 
nant  point  away  from  the  droplet  surface. 

The  large  surface  blowing  also  cause  the  flow  separa¬ 
tion  to  occur  early  since  the  vertical  interaction  of  mass 
flux  momentum  with  the  boundary-layer  flow  makes  the 
pressure  gradient  more  favorable  for  flow  separation. 
The  predicted  separation  angle  is  108®,  measured  from 
the  front  stagnation  point,  while  the  rigid  sphere  pre¬ 
dicts  129®  [21].  A  recirculation  wake,  detached  to  the 
surface,  is  clearly  shown.  The  wake  moves  backward 
in  downstream  direction  as  the  vaporization  increases. 
As  a  result,  the  increase  in  wake  length  will  be  also 
appreciable.  The  current  calculation  predicts  the  wake 
length  is  about  1.25  droplet  diameter.  For  the  same 
Reynolds  number,  the  sphere  without  surface  blowing 
has  a  wake  length  that  is  about  0.8  droplet  diameter 
[21]. 

We  have  learned  in  all  calculations  that  the  LOX 
droplet  surface  reaches  the  boiling  temperature,  which 
is  about  120  k  for  the  case  of  10  atms  cimbient  pres¬ 
sure,  very  quickly.  However,  the  droplet  core  remains 
cold  since  the  conductivity  is  low  and  at  early  time, 
the  circulation  strength  is  not  strong  enough  to  convect 
muchoe  energy.  All  the  available  heat  transfer  from  the 
gas  phase  is  used  to  heat  the  surface  and  also  utilized 
for  latent  heat  of  vaporization.  The  surface  boiling  per¬ 
sists  throughout  the  droplet  lifetime.  The  very  rich  and 
cold  oxygen  vapor  surrounding  the  droplet  is  depicted 
in  Figures  3  2ind  4.  The  large  gradient  occurs  at  the 
front  stagnant  region  as  we  expect. 

The  internal  circulation  of  the  droplet  is  evidenced 
in  Figure  5.  Since  a  stream  function  is  not  avail¬ 
able  in  this  unsteady,  compressible  problem,  we  have 
to  construct  the  “quasi-steady  stream  function'’,  which 
strictly  speaking  does  not  satisfy  the  continuity  equa,- 
tion  by  definition,  in  order  to  demonstrate  the  circula¬ 
tion.  Note  that  since  the  droplet  surface  moves  as  the 
vaporization  persists,  the  liquid  normal  velocity  at  sur¬ 
face  is  not  identically  zero.  Hence,  the  droplet  surface 
is  not  necessarily  a  streamline.  Due  to  the  large  viscos¬ 
ity  ratio,  the  time  needed  to  develop  a  spherical  vortex 
is  much  longer  than  for  the  hyrdrocarbon-fuel  droplet 
case. 

The  isotherms  of  droplet  at  two  different  times  are 
shown  in  Figure  6.  The  heat  conduction  mode  domi¬ 
nates  a  significant  portion  of  the  early  droplet  lifetime. 
The  convection  mode  then  takes  over  after  3  hydrody¬ 
namic  diffusion  times. 

The  pressure  distribution  along  azimuthal  direction 
is  displayed  in  Figure  7.  Due  to  the  momentum  dissi¬ 
pation  caused  by  the  surface  blowing,  the  pressure  re¬ 
covery  at  the  rear  portion  of  droplet  is  poor.  The  shear 
stress  distribution  is  illustrated  in  Figure  8.  The  large 


recirculation  wake  induced  by  the  early  separation  of 
the  flow  has  pushed  the  majcimum  shear  stress  point  to 
forward  of  the  equator.  The  asymmetric  distribution  of 
shear  stress  is  attributed  to  the  high  surface  blowing 

Figure  9  presents  the  surface  Nusselt  number  distri¬ 
bution.  By  a  pure  diffusion  analysis,  it  can  be  shown 
that  the  Nusselt  number  for  a  stagnant  sphere  is  2. 
The  results  indicate  that  the  Nusselt  numbers  in  the 
present  study  are  well  below  the  stagnant  value  even  the 
Reynolds  number  is  intermediate  high.  Similar  trends 
for  the  Sherwood  numbers  are  also  identified.  Due  to 
the  continuous  surface  boiling,  the  high  blowing  veloc¬ 
ity  has  significantly  impeded  the  heat  and  mass  trans¬ 
port  processes  of  the  LOX  droplet. 

The  normal  velocity  is  mainly  determined  by  the  sur¬ 
face  mass  fraction  and  mass-fraction  gradient  which  are 
associated  with  the  surface  temperature,  Nusselt  num¬ 
ber  cind  Sherwood  number  distributions.  Hence,  the 
distribution  of  normal  velocity  shown  in  Figure  10  re¬ 
sembles  those  of  surface  temperature  and  Nusselt  num¬ 
bers;  it  progressively  decreases  along  the  downstream 
direction,  then  slightly  increases  after  flow  separation. 
In  fact,  the  oxygen  vapor  at  the  recirculating  zone  is  rel¬ 
atively  lean  due  to  the  shape  effect  of  the  droplet.  As 
a  result,  the  blowing  velocity  increases  in  this  region. 

The  time  variations  of  total  drag  and  its  three  com¬ 
ponents,  pressure,  friction  and  thrust  drag,  are  illus¬ 
trated  in  Figure  11.  The  friction  drag  contributes  less 
than  20%  of  the  total  drag.  Also,  the  friction  drag  re¬ 
mains  almost  constant  throughout  our  calculation  since 
no  change  of  surface  condition  occurs  after  the  LOX 
droplet-surface  reaches  the  boiling  point.  The  slight  in¬ 
crease  in  pressure  drag  is  offset  by  the  increase  in  thrust 
drag;  as  a  result,  the  total  drag  almost  levels  off  even 
though  the  reduction  in  Reynolds  number  continues. 

The  transient  variation  of  Nusselt  number  is  demon¬ 
strated  in  Figure  12.  Care  must  be  taken  when  the 
convective  effect  is  taken  into  account  for  the  noncon- 
vective  Nusselt  or  Sherwood  numbers  in  a  simplified 
film-theory  model.  The  conventional  Frossling  corre¬ 
lation,  which  is  a  function  of  Reynolds  number  and 
Prandtl  (or  Schmidt  )  number,  may  yield  a  significant 
error.  Our  results  indicate  that  the  Nusselt  and  Sher¬ 
wood  numbers  are  much  less  than  their  stagnant  val- 
ues(2). 

Figure  13  presents  the  changes  of  droplet  heating 
rate,  averaged  surface  mass  fraction,  temperature,  and 
averaged  volumetric  temperature  (nondimensionalized 
by  the  boiling  temperature).  Due  to  the  persisting  sur¬ 
face  boiling,  the  energy  distributed  to  the  droplet  heat¬ 
ing  is  moderately  low.  Most  of  the  available  energy  for 
LOX  droplet  has  been  used  for  enthalpy  of  vaporiza¬ 
tion.  The  surface  temperature  peaks  very  quickly.  The 
volumetric  temperature  increases  slowly  with  time. 

A  case  with  —  20  atms  is  also  performed.  The 
qualitative  behavior  is  similar  to  that  of  case  of  = 
10.  Since  our  supercomputer  resource  is  very  limited. 
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we  decide  to  terminate  the  calculation  when  30%  of  the 
mass  has  vaporized.  Our  results  indicate  that  the  LOX 
droplet  vaporizes  much  faster  in  the  high  pressure  case 
as  shown  in  Figure  14.  The  diffusion  time  scales  in  the 
high  pressure  environment  are  shorter  than  in  the  low 
pressure  environment.  The  period  for  transient  droplet 
heating  and  development  of  the  internal  circulation  de¬ 
creases  as  pressure(  and  density)  increases  and  thus  the 
droplet  lifetime  reduces.  The  preliminary  results  also 
indicate  that  the  drag  coefficients  is  relatively  insensi¬ 
tive  to  the  variation  of  pressure.  However,  the  Nusselt 
and  Sherwood  numbers  increase  as  pressure  increases. 

We  have  experienced  numerical  difficulties  with  the 
high  pressure  cases  (p^  >  40  atms).  The  surface  tem¬ 
perature  shoots  above  the  critical  temperature  at  the 
very  ejirly  time  of  the  computation.  The  liquid  phase 
pressure  solver  fails  when  the  interface  solutions  be¬ 
come  disarraied  due  to  the  overshoot  of  surface  tem¬ 
perature.  In  our  first  analysis,  two  assumptions  of 
our  model  turn  out  to  be  critical  in  the  high  pressure 
case.  In  the  present  model,  there  is  an  inconsistency  of 
equation  of  state  employed  at  the  gas/liquid  interface. 
We  use  the  Redlich-Kwong  equation  with  mixing  rules 
to  determine  the  mole  fraction  and  then  the  density, 
while  in  the  gas-phase  pressure  equation,  we  employed 
the  ideal-gas  law  to  predict  the  density.  The  inconsis¬ 
tency  grows  as  the  pressure  increases.  The  second  bad 
assumption  is  that  the  solubility  of  fuel-vapor  in  the 
LOX  droplet  is  neglected.  Usually,  the  consideration  of 
real  gas  effects  will  predict  higher  mole  fraction  than 
Raoult’s  law  predicts.  When  the  droplet  approaches 
the  near  critical  state,  our  fugacity  and  mole  fraction  it¬ 
eration  procedure  predicts  a  mole  fraction  greater  than 
unity  and  destroys  the  balanced  conditions  at  inter¬ 
face.  The  omission  of  absorption  of  fuel-vapor  into 
LOX  droplet  also  drastically  affects  the  estimation  of 
enthalpy  of  evaporation. 

5  Concluding  Remarks  and  Fu¬ 
ture  Work 

We  have  conducted  a  very  detailed  study  of  a  LOX 
droplet  vaporizing  in  a  convective  environment.  The 
primitive- variable  approach  has  been  employed  to  solve 
the  complete  set  of  Navier-stokes,  energy,  and  species 
equations.  The  preliminary  results  for  low  pressure  case 
have  provided  us  a  clear  picture  of  the  transport  pro¬ 
cesses.  The  high  surface  blowing  velocity  has  a  signifi¬ 
cant  effect  to  the  flow  structure  and  modifies  the  flow 
separation  angle  and  wake  length  etc.  The  drag  coef¬ 
ficient,  Nusselt  and  Sherwood  numbers  are  reduced  to 
below  their  corresponding  stagnant  values.  The  sur¬ 
face  boiling  occurs  immediately  after  the  LOX  droplet 
is  injected  to  combustor  and  persists  throughout  the 
droplet  lifetime.  The  LOX  droplet  would  spend  a  large 
portion  of  energy  in  providing  the  heat  of  vaporization 


rather  than  heating  the  internal  fluid.  .\s  a  result,  the 
transient  droplet  heating  is  responsible  for  the  unsteady 
thermal  behavior  of  the  LOX  droplet. 

This  research  by  no  means  is  completed.  The  cur¬ 
rent  model  apparently  is  insufficient  to  handle  the  LOX 
droplet  at  the  high  pressure  environment.  The  ideal  gas 
assumption  is  questionable  when  the  pressure  is  near 
the  critical  condition  where  the  interaction  of  other 
species  is  important.  There  is  a  lack  of  a  computa¬ 
tionally  efficient  equation  of  state  which  is  applicable 
for  a  wide  range  of  pressure  in  two  dimensional  cal¬ 
culations.  The  Redlich-Kwong  equation  is  accurate 
but  is  too  time-consuming  to  incorporate  with  a  pres¬ 
sure  correction  equation.  The  new  method  to  linearize 
the  nonlinear  equation  of  state  is  currently  underway. 
A  very  restrictive  assumption  that  ignores  the  solu¬ 
bility  of  the  fuel-vapor  in  the  liquid  phase  should  be 
released  in  our  future  computation.  There  is  also  a 
need  to  review  more  advanced  numerical  techniques  at 
the  present  high  transfer  number  and  high  density-ratio 
two-phase  flow.  The  ADI  originated  methods  currently 
employed  converges  slowly  at  the  early  computation. 
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GAS-PHASE  VELOCITY  VECTORS  LIQUID-PHASE  STREAM  FUNCTION 


Time  =  5.00  ,  Instantaneous  Reynolds  Number  =  83.85  Time  =  5.00  ,  Instantaneous  Reynolds  Number  =  83.85 

Figure  2:  Instantaneous  velocity  of  gas  phase  at  time  Figure  5;  Stream  function  of  liquid  phase  at  time  =  5 

=  5. 


GAS-PHASE  ISOTHERMS 

Contour  tainvat  4.42E+01  Min;  1.19E+02Max:  l.OOE+03 


Time  =  5.00  ,  Instantaneous  Reynolds  Number  =  83.85 


LIQUID-PHASE  ISOTHERMS 

Coniour  Inlerval:  1.38E+00  Min:  1.00E+02Max;  1.21E+02 


-10  0.0  I.O 

Time  =  2.51  ,  Instantaneous  Reynolds  Number  —  91.99 


Figure  3:  Temperature  contour  plot  at  time  =  5. 


OXYGEN  MASS  FRACTION 


Contour  Iniervab  6.73E-02  Min:  -9.34E-03  Max;  l.OOE+00 


A. 9  0.0  4.9 

Time  =  5.00  ,  Instantaneous  Reynolds  Number  =  83.85 


Figure  4;  Mass  fraction  contour  plot  at  time  =  5. 


LIQUID-PHASE  ISOTHERMS 


Contour  Interval;  I.OIE+OO  Min:  I.06E+02  Max:  1.21E'MD2 


Time  =  5.00  ,  Instantaneous  Reynolds  .Number  =  83.85 

Figure  6:  Transient  history  of  LOX  droplet  heating  at 
two  different  times. 
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_  tkM  S  0.7,  H»(l)  =  07  local  min.  at  80.0 

_ llmas  2.6.  Ha(l)  =  91.7,localmln.al  66.8 

_ ttnos  Si.  Ra(l)  =  83.0,  local  min.  al  80.2 

_  tkno  =  7i.  Ra(t)  =  7S6,  local  min.  at  03.0 

_  tkno  =  10.2,  Ra<t)  =  67.4,  local  min.  at  94.1 


_  tlmo=  0.7, na<t)  =  97i,a4i. at  110.7,  local  min. all 3S.6 

_ tknas  2.0,  na(t)  =  01.7,  a4>.  at  lOOi,  local  min.  at  138.8 

_ llma=  5i,Ra<l)  =  834l,a.p.  at  lOOi,  local  min.  at  141.1 

_  tkno=  7i,Ra<t)  =  75i.a4i.  at  110.1,  local  min.  at  141i 

_  timas  10.2,  Ra(t)  =  67.4,  aji.  at  112.0,  local  min.  at  142i 


ANGULAR  POSITION  (DEGREE) 


0.  18.  36.  54.  72.  90.  108.  126.  144.  162.  180. 

ANGULAR  POSITION  (DEGREE) 


Fig.  7  Surface  pressure  distribution  at  different  times. 

_  tbiM  =  0.7,  Hi(t)  s  97  J,  a.p.  at  110.7,  atroao^  at  110.7 

_ tkoas  2.8,  Ra(tls  91.7, 841.  at  10041,  atratteO  at  10841 

_ tkiw:  Si,  Ho(t)  s  83.0, 8.p.  at  109i,  atioaosO  at  109i 

_ tkiw  =  7J,  Bo(t)  r  7Si,  84).  at  110.1,  MraaacO  at  110.1 

_  timos  10.Z  B8<t)  =  87.4, 8.p.  at  112.0, 8ti888=0  at  11241 


0.  18.  36.  54.  72.  90.  108.  126.  144.  I6Z  ISO- 

ANGULAR  POSITION  (DEGREE) 


Fig.  9  Surface  Nusselt  number  distribution  at  different  times. 

_  tkiMc  0.7,  Ra(t)s97i,  84).  at  llOi,  local  min.  8113041 

_ ttaws  2.6,  Haft)  s  01.7, 8.p.  at  lOoi,  local  min.  at  135.1 

_ tbiws  Si,  Baft)  s  8341, 84>.  at  lOOi,  local  min.  at  1368 

_  tfenos  78,  Baft)  s7S8,8.p.  at  110.1,  local  min.  at  141.6 

_  tfcnoc  lOi,  Haft)  =  67.4,  a.p.  at  11Z0,  local  min.  at  1428 


0.  18.  36.  54,  72.  90.  108.  126.  144.  I6Z  180. 

ANGULAR  POSITION  (DEGREE) 


Fig.  8  Surface  shear  stress  distribution  at  different  times. 


Fig.  10  Surface  normal  velocity  distribution  at  different  timei. 


1.  Tolsl  Drag,  p  «  10  aims. 

2.  Prtssura  OniO  Componant. 


0.0  1.0  2.1  3.1  4.1  5.2  6.2  7.2  8.3  9.3  1( 


Gas  Hydrodynamic  Diffusion  Time  Scale 


-  1.  Oroplat  Hasting  Rata  (Ql/Og) 


0.0  1.0  2.1  3.1  4.1  5.2  6.2  11  83  9.3  10.3 


Gas  Hydrodynamic  Diffusion  Time  Scale 


Figure  11:  Time  variations  of  total  drag  and  three  drag  Figure  13:  Time  variations  of  droplet  heating  rate,  av- 

components.  eraged  surface  mass  fraction,  temperature,  and  aver¬ 

aged  volumetric  temperature. 


0.0  1.0  2.1  3.1  4.1  5.2  6,2  7.2  8.3  9.3  10.3 


Gas  Hydrodynamic  Diffusion  Time  Scale 


-  1.  TgslOOO  K,  TIslOO  K.  LOX,  R»*100,  p.10  Mm. 

- 2.  70=1000  K.  71=100  K,  LOX,  RtelOO,  p=20  »tm. 

3.  Mass,  p«l0  aim 

4.  Mass,  pe20  atm 


0.(X)  1.03  2.07  3.10  4.M  5,17  021  7  24  8,28  9,31  10.35 

GAS  HYDRODYNAMIC  DIFFUSION  TIME  SCALE 


igure  12:  Time  variations  of  averaged  surface  Nusselt 
umber. 


Figure  14:  Time  variations  of  radius  and  mass  for  the 
cases  of  different  ambient  pressures. 
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